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ABSTRACT 


Control of the motions and vibrations of large space structures require the know- 
ledge of state values that mav not be available due either to inability to measure the 
states or, the high cost of the sensors to measure the required states. One solution is the 
use of an observer to estimate the states from limited sensor input. 

The physical characteristics of large space structures and the enviroment they oper- 
ate in Will cause large amounts of noise in the measurements. The obvious observer for 
such an enviroment is the Kalman Filter which is specifically designed to produce opti- 
mal estimates in a noisvV enviroment. 

A straightforward application of the Kalman Filter will be examined utilizing a 
steady state Kalman gain matrix. The observer performance will be examined in both 


matched filter’plant and reduced order filter configurations. 


THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While every effort has been made, within 
the time available, to ensure that the programs are free of computational and logic er- 
rors, they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION 
A. BACKGROUND 


The advent of large space structures poses a number of problems for the control 
engineer. Previously, the objects put into space could be treated as rigid bodies so that 
a single three axis sensor package could be used to tell the motion of all components. 
The large space structures will not be rigid, instead they will have considerable flexibility 
and multiple modes of vibration [Ref. 1: p. 51] 

Control of the structure’s attitude and vibrations requires knowing the motions of 
the components. One approach would be to heavily instrument the space structure, but 
weight and cost make this approach impractical. An alternative is to use a limited 
number of sensors to measure only certain states and to deduce the other required states 
by use of an observer algorithm. 

This thesis will address the production of estimates of the states needed for control 
of the structure. The model used will be a early design study by McDonnell Douglas 
Astronautics for a dual keel space station. The techniques and problems of observation 


for this model are generic to all large space structures. 


B. PROBLEM STATEMENT 

Design of an observer for estimating the states of a large space structure breaks 
down into several steps. First, a matematical model is developed for the system behavior 
over time. Modal analysis is used to form a system composed of decoupled second order 
differential equations. The use of decoupled equations allows a reduced order model to 
be generated by truncating the number of modal equations. A reduced order model will 
have all of the same mathematical qualities (and problems) but reduces the amount of 
time and computer resources required to do simulation. 

Second, the observer is designed. The observer 1s designed to obtain a minimum 
variance estimate of the desired state values from the measurements. 

Third, the observer is simulated to verify performance. Simulation runs of both a 
matched observer/plant system and a reduced order observer are employed. That is, the 
system is run where the observer is used to estimate all of the plant states and run where 


there are more plant states than the observer estimates. 


Fourth, results are analysed and conclusions drawn based on these results. 
Reccomendations for further areas of research are suggested based on the results and 


conclusions. 


C. ORGANIZATION 

The model of the space station is developed in Chapter II. The modal model was 
developed using modal analysis and discretized to form the discrete-time state equations. 
The data for this model was from an early design study by McDonnell Douglas 
Astronautics Company for a dual-keel space station. The observer and its equations are 
developed in Chapter III. Chapter IV is the simulation runs of the observer versus the 


plant. Chapter V presents conculsions and recommendations for further research. 


Il. MATHEMATICAL MODEL 
A. INTRODUCTION 


Prior to the proposed space station almost all of the objects put into space could 
be treated as simple rigid bodies for the purpose of mathematical modelling of their 
motions. The design constraints imposed by the high cost of lifting mass to orbit dic- 
tates a light, open structure with considerable flexing. Large space structures such as the 
space station, therefore, cannot be treated as rigid bodies. The structure is in fact lightly 
damped with multiple natural frequencies. The result is a structure that will vibrate for 
considerable periods of time whenever external forces are applied. 

The space station structure can be modeled as an n-DOF (degree of freedom) svstem 
consisting of n masses, springs, and dashpots [Ref. 2: p. 173-176]. This straight forward 
modelling of the coupled masses produces a system of unworkable complexity. As a 
result, the system will be modelled in terms of the structures natural modes of vibration. 
The resulting system, while still complex, is at least workable. 

The model will be developed in two steps. The first will be to generate the 
continous-time model of the natural modes. The second will yield the discrete-time 


model, developed from the first model, for use in the simulation. 


B. MODAL MODEL 

The space station structure can be modeled as a system of discrete masses coupled 
by springs and dashpots. The major mechanism of damping in the structure 1s structural 
damping, the internal dissipation of energy within the members, as the structure vibrates. 
Structural damping can be shown to be equivalent to viscous damping and this equiv- 
alency is used in the model [Ref. 2: p. 72-73]. 


The energy dissipated by structural damping 1s: 


Wd = aX? (1) 
Wd = energy dissipated by structural damping 
o = constant (force/displacement) 
X = displacement 


The energy dissipated by viscous damping 1s: 


Wy = ncwX? (2) 


We can equate the two 
mCegtoX° = aX? (3) 


yielding an equivalent viscous damping coefficient: 


Ca, = — (4) 
The second order differential equation for a single viscously damped mass 1s: 


mx + cx + kx = F(t) (S) 


Substituting C,, for c 
mx + x + kx = Fs) (6) 


For multiple mass systems C,, becomes oa where w, is the natural frequency of vi- 
bration. 

The displacement of masses can be represented by the second order matrix differ- 
ential equation [Ref. 3: p. 3-9], 


a ees 
Mgt) + pp Kali) + Kgl) = F (2) (7) 

g = coordinate vector 
M = system mass matrix (diagonal) 
&K = equivalent damping 
d = damping coefficient 
w, = frequency of oscilation of the system 
K = svmmetric system stiffness matrix 
F(t) = system forcing function 


The above equation represents a system of second order differential equations cou- 
pled through the stiffness matrix. Decoupling can be done by expressing q in terms of 
natural modes of vibration. The process is called modal analysis. The independent dif- 
ferential equations can then be treated individually. The modal equations are derived 
below. 


First, the undamped, homogeneous form of Eq. (7) 


M4(2) + Kg(r) = 0 (8) 


is solved. Let 


q(t) = Ax sin(wr + ©) (9) 
q(t) = Axw cos(wi+ O) (10) 
q(t) = — Axw’ sin(wr + ©) (11) 


Stosmiuune bq. (9) and Eq. (10) into Eq. (11) 
[ — w°M + KJAx sin(wr + ©) = 0 (12) 
This equation has a non-trivial solution for all time if and only iff 
[K — wo M]x =0 (13) 


Equation (12) has n combinations of x (natural mode shapes) and @ (natural fre- 


quencies) as solutions. These can be grouped into matrices: 


X = [x)x5...20,)7 (14) 
Q? = diagl a, 97--Oon] (15) 

which satisfy the equation: 
KX = Q°MX (16) 


Several useful relations can be derived from Eq. (16). Premultiplying Eq. (16) by 
ee 


X’KX = 2°X’MX (17) 


The eigenvectors can be normalized 


X’MX =I (18) 
which yields 
X’KX =’ (19) 


The equations of motion can be uncoupled through the linear transformation of the 


coordinate system 


nr 


g(t) = > x(t) = Xn(0) (20) 


i=t 
X = modal matrix 
n = maximum number of degrees of freedom 


y(t) = transformed coordinate vector 


Application of the transformation to the system Eq. (7) yields 
X™MXij(1) + or X7KXi (0) + X7KXn(1) = X7F() (21) 


Using Eq.(18) and Eq. (19) 


X/MXij(2) = li(d) = 7 (22) 
A eTus- Hao: {i 
o, X’KX7(0) = ain 741) = dQ (23) 
X7KXn(1) = Q°7 (24) 
therefore 
7 + dQy + Q°n =X'F (25) 


Equation (25) is the modal model of uncoupled second order differential equations. The 


motion of the structure can be found from the modal amplitudes, y(z), using Eq. (20). 


C. DISCRETE-TIME MODEL 
The discrete-time state space model is found by solving the continuous-time 


equations. The ith equation of motion is 


Hit) + dedoit(t) + woin(t) = Xj F(0) (26) 
X? = transpose of the ith mode shape vector 
F(t) = torquing force applied at a point 


The homogeneous solution (F(r) = 0) for Eq. (26) is (Ref. 4: p. 475-476] 


ndt) = Cie” cos(wgt) + Cpe” sin(w,t) (27) 


where 


= > (28) 


OTF 05, = y (29) 
The constants in Eq. (27) can be found by taking the derivative 
n(t) =(Cywg— Cyy)e"' cos(wyt) — (Ci\aq— Coy)e” sin(w,) (30) 


and evaluating at r= 0 


(0) = C, 2 (31) 
n(O) = C,wg— Cyy 2) 

Solving for C, and C, 
C, = 7(0) (33) 


#0) | nO 


Wa Wa 


e 1 0 I[n(0) 
=| , 3 
| [t to ae 


Rewriting Eq. (27) and Eq. (30) in matrix form 


n(t) Z e”' cos(w,t) e”’ sin(w gt) A 36 
n(t)| |e Ly cos(wgt) + wg sin(w,t)] e7”[w, cos(w,t) — y sin(w,t)] oe 


Substituting Eq. (35) into Eq. (36), the solution can be written in terms of the initial 


C, — 


In matrix form 


conditions 





PA _ eM cos(@al) + —- sin(w ft) | a poe sin(@¢!) on (37) 


1(t) ~ - on d~”" sin(w gt) eT cos(wal) — a sin(wgt)] {| 7(0) 
Letting 
t 
X() = iw i (38) 
nt) 
and 
eee l =e 
we e_y{ cos(aggt) + Ge sin(wge)] aa ms mee 39) 
. — ion e~”* sin(w tf) eT cos(w,t) — a sin(wyt)] 
the solution can be written as 
X(t) = A,(t)X;(0) (40) 


where A, is the state transition matrix of the ith mode. The non-homogeneous solution 


1S 
X{) = A{1)X{0) + B.x/ F(0) (41) 


where the discrete-time input matrix, for constant F, is given by 


B= | Bore (42) 
0 


and T =[0 1]? is the input matrix for the continuous-time system, and T is the sampling 
time. Solving Eq. (42) yields 


] 
B;= @ 


y 





[l1—e ” cos(wyt) — = e7 ’ sin(w 4t)] 


; (43) 


=v ieee 
map sin(@ 4) 


The discrete-time state equation for the ith equation of motion can be written as 
X(AT + 1) = A{T)X{AT) + B(T)x/ F(AT) (44) 


where A, and B, are evaluated atr=T7. Here, 


X, = vector of the ith modal amplitude and the ith modal velocity 


A, = ith state transition matrix 


B, = 1th input vector 
x7 = transpose of the ith mode shape vector 
F = distrubuted force on the plant 


= sampling time 


time index 


~ 
I 


Equation (44) can be expanded to include the disturbance input, w({kT) : 
X(AT + 1) = A(T)X{KT) + B(T)x; [F(AT) + (kT) (45) 


Equation (45) is the discrete-time mathematical model describing the motion of the 


structure in terms of its natural modes of vibration. 


Il. THE OBSERVER 
A. INTRODUCTION 


The observer design will be required to estimate the modal states in a noisy 
enviroment. Kalman filtering is the most widely used technique for accomplishing the 
production of state estimates in a noisy enviroment [Ref. 5: p.159]. The steady state 
Kalman filter was selected to minimize the computations during the actual plant ob- 
Server operation. Use of a steady-state gain matrix for the observer allows the matrix 
to be computed separately from the operational observer, reducing the computer power 


required for the observer and allowing the algorithm to operate more rapidly. 


B. KALMAN FILTER EQUATIONS 
The discrete Kalman filter provides state estimates for the following dynamic sytem 
[Ref. 5: p. 159-162], 


X(A + 1) = AN(A) + BU(A) + BnW(k) (46) 
Y¥(A + 1) = CX(A + 1) 4+ V(K +1) (47) 
X =nx1I1_-_ state vector 
Ue — Px) control vecton 
Woe= 1 xX lPplanmneisexccron 
Y = mx1_= measurement vector 
VY =mx1_= measurement noise vector 
A =nXn._ State transition matrix 
B =nxp control input matrix 
Bn = nxr_ plant noise input matrix 
CC) = inn Smeastrementmannn 


The plant noise vector W(k) is gaussian white noise with 
E{W(k)} =0 (48) 


E{(W(A)W’(4)} = Q (49) 


10 


for all A =0,1,2,..., and Q is a positive semi-definite r x r matrix. V(k) is gaussian white 
noise with 


E{V(k)} =0 (50) 
E{V(A)V(kK)} =R (51) 


for all A =0,1,2,... , and R is a positive definite m x m matrix. The two random processes 
W(k) and V(k) are assumed to be independent, so that 


E{VQ)W(A)} = 0 (52) 


for all j= 1,2,..., and A =0,1,2,... The intial state X(0) is assumed to be a gaussian ran- 


dom vector with 
E{X(0)} = 0 (53) 


It is assumed that X(0) 1s independent of W(A) and V(&). 
The optimal estimate of X(k + 1) 1s denoted X(k +1|/k+1). The Kalman filter is 


designed to minimize 
J=E{(X(k + 1)- X(k +1]k+ 1 [X(k +1)- X(k +1/k+1)}} (54) 


The recursive realtions for generating X(k +1{k+1) are 


K(k 41] kb) = AX(k | &) + BU(K) (55) 
K(k + IVR + 1) = XK 411 R$ GR + IVR + I — CXR +1 0) (56) 


fore~ — 0,1,2,... , where X(0 | 0) = 0. X(0 | 0) is set equal to zero since the expectation of 
X(0) is Zero. 
G(k + 1) is ann X m matrix , called the Kalman Gain Matrix which 1s specified by 


the realtions: 


P(A +1] k) =AP(K| AA’ + BnQ(A)Bn’ (57) 
Gk +1 =P(AF1/ AC ICP(K +1) 4074+ R(K4+ D7! (58) 
P(kK+ LE; A+ 1) =[I— GK + ICIP(A + 1 4) (59) 


P(x | A) is the covariance matrix of the error between the states and their estimates 

P(k | k) = E{EX(k) — X(k | AEX) — X(N) (60) 
Since We are using the steady state gains the choice of P(0| 0) is irrelevant. P(0|0) is 
initialized to zero in the gain derivation program for simplicity [Ref. 6: p. 139-140]. 


C. STEADY-STATE SOLUTION 


If Equations (57), (58), and (59) are repeatedly iterated, G(k + 1) will converge to a 
steady state value |Ref. 7: p, 263): 


G,,= lim G(k+1) (61) 


k > c”0 


The values of G,, (or G) can be substituted into Eq. (56) making the steady state 
Kalman filter 


X(k +1] =AX(k| + BU(K) (62) 
Kk+ 1 A+ I =X(K 41/4) 4 G6(Y(K 4 I) — CRA 41/1] (63) 


D. OBSERVER PERFORMANCE 
The performance of an observer is judged by how accurately and rapidly it estimates 
the desired states. The performance measure of the observer as a whole is shown in 


equation (54). The normalized performance of the observer for individual states is 
J, = El(x) — 5) VEL] (64) 


which can be found using Eq. (65) 


J;= » ett — (k)) T, + » F007, (65) 
k=0 k=0 
a = performance measure for the ith state 
x(k) = value of the irh state at k 
x(k) = observer estimatecof ith state at k 


ie = sample interval 


A normalized performance measure is used to aid comparison of the performance of the 
observer in estimating various states. From Eq. (65) it can be shown that if x(k) = 0 for 
all k =0,1,2,... that J; would be unity. Therefore, the better the performance of the ob- 
server, the smaller the fraction of one J, will be. 


13 


IV. SIMULATION 


A. INTRODUCTION 


The objectives of the simulation were to 


e determine the sensitivity of the observer performance and settling time to changes 
in the ratio of plant noise to measurement noise, 


e determine the effect on observer performance and settling time of increasing the 
number of modes observed in the matched plant/observer, and 


e determine the performance for the reduced order observer. 


B. PLANT AND OBSERVER DATA 

The dynamic model is a truncated form of a preliminary space station configuration; 
the phase II dual keel structure.! The full model consists of an infinite number of natural 
modes but this was restricted to the first ten active modes for this study due to limita- 
tions On computer resources. As will be shown reasonable data can be obtained with this 


simplification in examining the observer performance. 


C. SIMULATION PROGRAMS 

The simulation was broken down into two segments due to the large memory and 
computational time requirements. The first program computed the steady state observer 
gain matrix (G). The second program ran the observer and the plant when the plant 
was subjected to an impulse excitation. 

The steady state observer gain matrix (G) was obtained by repeated iteration of 
equations (57), (58), and (59). The equations were were run until the values of the ma- 
trix changed by less than a set fraction. The following formula was used to check the 


changes in the gain matrix elements 
Agi; = Leifk + 1) — 8; fA) Sek + 1) (66) 


The program was terminated when og,, was less than 10-". 

The settling time for the estimates of the states to be within 2% of the actual states 
was determined by finding the eigenvalues of A - G*C then computing as follows [Ref. 
6: p. 139-143] 


1 The model for preliminary station configuration was provided courtesey of McDonnel 
Douglas Astronautics Company, 5301 Bolsa Avenue, Huntington Beach, CA 92647. 
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Ts = log(.02)/ log acc) (67) 


The expected error in the sensor, i.e., the standard deviation of the noise in the 
measurement, was choosen as 10° feet per sec. per sec. based on the natural frequencies 
in the structure and reasonable sensor sensitivity [Ref. 2: p. 79-80]. The expected plant 
noise was varied to find the range of ratios between plant and measurement noise that 
the filter would be effective. This approach was taken since the plant noise contributors 
are not currently well defined. 

The second program subjected the plant as modeled in Eq. (45) to an impulse 
excitation and then had the observer estimate the selected states using observer 
equations (62) and (63). Observer performance was computed using Eq. (65). 

A third program was used to find the contribution of unobserved modes to the noise 
in the Kalman observer. The program ran the plant subject to an impulse excitation and 
computed the product of the measurement matrix C times the unobserved modes of the 
State vector X(A) for a measure of the noise contributed bv the unobserved modes. 


The three programs are listed in the appendices. 


D. EFFECT OF PLANT TO MEASUREMENT NOISE RATIO ON OBSERVER 
PERFORMANCE 

The ratio of the variance of the plant noise (PN) to the variance of the measurement 
noise (MN) was found to have a strong effect on the Kalman Observer performance (J) 
and settling time (7,). Figures 1 through 6 show the observer performance for a 3 mode 
matched plant and observer system for progressively smaller PN/MN ratios.2 Figure 7 
shows the performance for the seventh mode (position) versus several values of PN’ MN. 
Figure § is the settling time versus the same PN/MN ratios. 

The figures show that, for all of the plotted performance values, the observer per- 
formance is at least marginally acceptable regardless of the PN/MN ratio. Decreasing 
the PN’/MWN ratio leads to an even more rapid degradation in observer performance. 


The settling times also rapidly increase as the PN/MN ratio decreases. 


2 Figures 1-6 and 9-15 show the performance measure for each mode. The bar for the mode 
position is immeadiately to the right of the numbered tick mark on the x-axis scale, the mode ve- 
locity is next to it immeadiately adjacent to the tick mark without a number. 


E. EFFECTS OF INCREASED MODES ON OBSERVER PERFORMANCE 

The matched plant/observer was run with increasing numbers of modes to see if 
there was an effect on observer performance (J) or settling time 7,. Figures 7 through 
17 are of observer performance for systems with increasing numbers of modes in the 
system being observed. Figure 18 is of settling time versus the number of modes in the 
system. The ratio of PN/MN was kept constant at PN/MN = 2.5 x 10°. 

The increasing of the number of modes for the matched plant/observer had neglible 
effect on the performance for the individual modes. The performance value for the 
modes was effectively constant. Settling times for the observers increased as the number 


of modes was increased. 


F. REDUCED ORDER KALMAN OBSERVER 

The Kalman Observer has been shown to be effective where the number of modes 
observed matches the number of modes in the plant. The Kalman Observer was then 
run with the one less mode observed than the number of modes in the plant. The gain 
matrix (G) from the matched system was used. The observer failed with the state esti- 
mates produced by the observer becoming excessively large and having settling times of 
hours vice minutes. Since the purpose of the observer was to provide estimates for use 
in controlling the plant the time delay makes the estimates unusable. 

The cause of the observer failure is apparent when you look at the last portion of 


Eq. (56) of the Kalman Observer 
G[Y(k + 1) — CX(k +1] 4)] (68) 


This portion of the observer equation is the correction of X(k+ 1|k) to produce 
Xk +1|k+1). The design of the Kalman observer is to produce an estimate despite 
the measurement noise but, with the reduced order filter there is additional unanticipated 
noise which causes over correction of the values of X leading to the state estimates being 
excessively large and settling times being too long. This can be shown by examining 


What composes Y(k + 1) — CX(k + 1| 4) 


16 


C]| --- |—Cl x,(k) (69) 


C times the state x,_,(k) through x,(A) is unanticipated noise so if 


x,(4) x,(A) 
x(k) x,(k) 
1 1 
c —C =) 70 
| | mo) 
aes (K) Be ACs) 
Xm(k) Xm(K) 


the remaining portion of the C matrix times the modal states is an equivalent noise. 
Table (1) shows the growth of the unanticipated noise in the filter as the number 
of unobserved modes in the plant grows. Table (2) shows the individual contributions 
of the individual modes when left unobserved. Table (1) shows that the unanticipated 
noise is much larger than that expected by the filter (10-*). Table (2) shows that there 
are modes that do not markedly contribute to the noise and that they might successfully 
be left unobserved if the measurement noise estimate was already much larger than these 


noise sources. 


Table 1.5 CUMULATIVE UNANTICIPATED NOISE FROM UNOBSERVED 
MODES 


ees aas [6208 aon 
[soar _arse7ag (|e | —092.50 
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Table 2, UNANTICIPATED NOISE FROM UNOBSERVED MODES BY 
MODE 
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V. CONCLUSIONS AND RECOMMENDATIONS 
A. CONCLUSIONS 


Simulations runs showed that a matched plant/observer can work if the following 
criterions are meet: 


¢ The ratio of plant noise to measurement noise is sufficiently high to produce a us- 
able settling time. 


e Sufficient computational power is available to run the matched observer. The 
amount of memory and number of computation goes up as the number of modes 
observed increases. 

Utilizing a reduce order observer for an arbitrarily selected set of modes is not fea- 
sible. The non-observed modes add so much noise to the system that settling times and 
observer performance are so poor as to render the observer useless for obtaining state 


values for plant control. 


B. RECOMENDATIONS 
The work on the Kalman Observer for Large Space Structures lead to the following 


recommendations for further research: 


e Identify those modes that contribute the largest noise to the Kalman observer and 
set the observer to estimate these states in addition to those required for plant 
control. A possible method for identifving the modes that contibute the largest 
noise to the observer would be the Karhunen-Loeve expansion. 


¢ Modifving the plant/observer to model the use of sensors at additional positions 
to see if the increase in the data rate will help decrease settling time. 


¢ Modify the model to incorporate noise injection into more than one location. The 
current model has noise injected at only one position, a useful simplification for 
initial analysis but not realistic. 
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APPENDIX A. KALMAN GAIN MATRIX GENERATION PROGRAM 


sSevevedededesedetevetetetedosetetete sede tedesete fete Fehe Soh PoE RICCI RTA RE RAR RAR AERTS 
oes cae GGAIN Sesvevedveve 
weet ADAPTED TO RUN KALMAN FILTER AND COMPUTE THE = ¥iiier 
we = G MATRIX BY ITERATION STOPPING WHEN THE Feeseseve 


week = =THE MATRIX GOES TO STEADY STATE Feveverese 
Sevedesevesestesodedevedededededectetetedetototetcde tet dedcsedede kde cede Feu RRR RRRRARERRRREERAR RR 


Sevedevevedestevovessdetedededededetedstedesesete de tesecededede Ref RRR RR RRA RRRRERRIER 


aussi VARIABLE DECLARATIONS KIKI 
dered dedede vere deve dete deter de FePeNete RII RI ILLICIT RII IARI AISI 


EXTERNAL STMTRX,DLINRG,EXCMS, DEVCRG 

CHARACTER*6 NAM 

CHARACTER*1 AGAIN, CORECT,RAGAIN 

INTEGER ROWN1, ROWN2 ,ROWN3, COUNT , NODE, MODE ,KQ,EMODE , SMODE ,R2M,C2M 
INTEGER CT,CF,KADJ,CFADJ, LOOP,PRNT,JJ,JK,N1,JR,KR,MR, ISEED ,M2 
INTEGER JL,J1,JM 

REAL LAMA(100), UGVEX(684, 100) ,RNODE,RMODE ,MIN 

REAL*8 PHI(2,2,100),GAMMA(2,100) ,EGT,GMA,WN,W1,X1T,X2T, TIME 
REAL*8 A(200,200) ,B(200,3),F(3, 50), IMPLSE,ENERGY 

REAL*8 COSW1T,SINW1T,X1(100) ,X2(100) ,COST( 100) 

REAL*8 DAMP, SAMPT, PI ,SAMPTM,SUM1,SUM2 ,SUM3,SUMC 

REAL*8 C(6,200), IDENT( 50, 50), RMN(6,6), QPN(3,3) 

REAL*8 PK( 50, 50), Y(6), BN(200,3) 

REAL*8 PNVARX, PNVARY, PNVARZ 

REAL*8 MNVX1, MNVY1, MNVZ1, SUM, BQBT(50,50) 

REAL*8 TMP1( 50,3), TMP2(3,3), TMP3( 50, 50) 

REAL*8 PK1( 50, 50),G( 50,3) 

REAL*8 DY(3), ES,ED,ESUM,CGN, PRT 

REAL*8 SF, N9, TCHK, ACHK, H1, H2, H3, H4, HS, H6 

REAL*8 AGC(100, 100) 


COMPLEX*8 EVAL( 100), EVEC (100,100) 


tMieledskdkdchkhcchicsiiccnKkR Raa aR aikckRRRaRR«,R ate 


Aededeaeae VARIABLE DEFINITIONS Hedededeye 
VedstisthhKehkKkenKKcKRKKKKececccmcdckcdkdcckikkeeRKKk 


STMTRX = SUBROUTINE EXTABLISHES STATE TRANSITION MATRICIES 
LAMA = VECTOR OF THE SQUARE OF THE NATURAL FREQUENCIES 
UGVEX = MODE POSITONS AND SLOPES OF THE NODAL POINTS 

PHI = STATE TRANSITION MATRICIES FOR EACH MODE 

GAMMA = INPUT TRANSITION MATRIX 

A = DIAGONAL MATRIX CONSISTING OF PHI 

B = INPUT MATRIX OF GAMMA AND CONTROL SLOPES 

DAMP = DAMPING FACTOR 

SAMPT = SAMPLING TIME 
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GMA00010 
GMA00020 
GMA00030 
GMA00040 
GMA00050 
GMA00060 
GMA00070 
GMA00080 
GMA00090 
GMA00100 
GMA00110 
GMA00120 
GMA00130 
GMA00140 
GMA00150 
GMA00160 
GMA00170 
GMA00180 
GMA00190 
GMA00200 
GMA00210 
GMA00220 
GMA00230 
GMA00240 
GMA00250 
GMA00260 
GMA00270 
GMA00280 
GMA00290 
GMA00300 
GMA00310 
GMA00320 
GMA00330 
GMA00340 
GMA00350 
GMA00360 
GMA00370 
GMA00380 
GMA00390 
GMA00400 
GMA00410 
GMA00420 
GMA00430 
GMA00440 
GMA00450 
GMA00460 
GMA00470 
GMA00480 
GMA00490 
GMA00500 
GMA00510 


Get C2 C2 CD CO C2 C7 Cd O99 9 9 9 a) a rea oe OO Oe C2 CP OC) CY Cy Ce CO Cr OC) OG) a es eG eral) G) G) GG) @ @ 


TCX, TCY, TCZ = CONTROL TORQUE VALUES 

ENERGY = TOTAL SYSTEM ENERGY 

IMPLSE = IMPULSE INPUT FUNCTION 

MIN = NUMBER OF MINUTES SYSTEM WILL BE OBSERVED 
SMODE = NUMBER OF STARTING MODE (INT) 

MODE = NUMBER OF MODES (INT) 

EMODE = NUMBER OF THE LAST MODE (INT) 

NODE = NUMBER OF THE NOISE INPUT MODE (INT) 
eve NOISE SLOPE LOCATIONS IN DATA MATRIX #% 
ROWN1 = X-SLOPE LOCATION 

ROWN2 = Y-SLOPE LOCATION 

ROWN3 = Z-SLOPE LOCATION 

C = OUTPUT MATRIX FOR Y 

IDENT = IDENTITY MATRIX 

RMN = MEASUREMENT NOISE COVARIANCE MATRIX 

QPN = PLANT NOISE COVARIANCE MATRIX 

PNVARX = PLANT NOISE X-SLOPE VARIANCE 


PNVARY = PLANT NOISE Y-SLOPE VARIANCE 
PNVARZ = PLANT NOISE Z-SLOPE VARIANCE 
MNVARX = MEASUREMENT NOISE X-SLOPE VARIANCE 
MNVARY = MEASUREMENT NOISE Y-SLOPE VARIANCE 
MNVARZ = MEASUREMENT NOISE Z-SLOPE VARIANCE 


ISEED = INITIALIZATION FOR RANDOM NUMBER GENERATOR 

XKAL = X MATRIX 

Y = OUTPUT MATRIX 

RNDM = RANDOM NUMBERS USED FOR WHITE NOISE IN MEASUREMENTS AND 
IN PLANT FORCES 

BN = B MATRIX TO MULTIPLY NOISE DISTURBANCES 

TNA, TINY ,INZ= NOISE TORQUES X,Y,Z SLOPES 

M2=2*MODE 


f, a = << aea/’aseselUlUlUlUlUlUlUlUlUlUlUUCc ‘Q(T A AA TPT LL af NEC CS TRA FF SL CF hb Lf tO OT UDWTClUCUWLhhltt ee ae eta awWawWawawa = = see = = = 
Fededodedetevevetsd tevededetetedotodendedededetstetede 


SAMPLE OF SPACE EXEC FILE 


THIS FILE MUST BEGIN IN COLUMN 1 AND RUN WITH THE FOLLOWING 
SEQUENCE FOR THE INITIAL RUN OF THE PROGRAM: 


FORTVS SPACE 
SPACE 
LOAD SPACE (START 


(COMPILES PROGRAM) 
CEXECUIES hE sG=r 1 LE) 
(LOADS AND EXECUTES PROGRAM) 


SUBSEQUENT PROGRAM RUNS CAN ELIMINATE "FORTVS SPACE" IF NO 
CHANGES HAVE BEEN MADE TO THE PROGRAM, AND CAN ELIMINATE 
RUNNING THE EXEC FILE. 


Pie tsk CHESTS [INPUT B (PERM 

Piecworsk UlILITY DATA (RECFM VS BLOCK 133°PERi 

FI 11 DISK CNTRL OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 
FI 13 DISK GAMMA OUTPUT (RECFM VS BLOCK 133 PERM 

FI 14 DISK MODE OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 
FI 16 DISK COST OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 
FI 17 DISK PRY OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 
FI 18 DISK ERROR DATA (RECFM F BLOCK 80 LRECL 80 PERM 
FI 19 DISK END FILE (RECFM F BLOCK 80 LRECL 80 PERM 

FI 20 DISK GMAT FILE (RECFM F BLOCK 80 LRECL 80 PERM 
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GMA00520 
GMA00530 
GMA00540 
GMA00550 
GMA00560 
GMA00570 
GMA00580 
GMA00590 
GMA00600 
GMA00610 
GMA00620 
GMA00630 
GMA00640 
GMA00650 
GMA00660 
GMA00670 
GMA00680 
GMA00690 
GMA00700 
GMA00710 
GMA00720 
MA00730 
GMA00740 
GMA00750 
GMA00760 
GMA00770 
GMA00780 
GMA00790 
GMA00800 
GMA00810 
GMA00820 
GMNA008 30 
GMA00840 
GMA00850 
GMA00860 
GMA00870 
GMA00880 
GMA00890 
GMA00900 
GMA00910 
GMA00920 
GMA00930 
GMA00940 
GMA00950 
GMA00960 
GMA00970 
GMA00980 
GMA00990 
GMA01000 
GMA01010 
GMA01020 
GMA01030 
GMA01040 
GMA01050 
GMA01060 
GMA01070 


aoanaan C2 €2;-0) 


aqaaqdan 


Lede te te teva dete te Tee Rede Be teRe VETERE EDEL EEL EEE ETE BEEREPEBEEEBEEBTEBREERBEEBRERRERE 
PARAMETER (JR=5243, KR=5397, MR=262139) 
MIN =1200.0 
WT=1. ODOO 
PI = 4.0D0 * ATAN(1. ODO) 
SesevevevevedeevevedetedesededesevedesekckdtedteteesetesKscdescklevededevekkhReRRdeseletccdtocksdkichk sede 
Yererererc READ LAMA AND UGVEX MATRICIES aedorevere 
sleslegleciocte cle sloslevoclovlovesiovioveseslesloveseveslesesevedesevevesesesevoveseselesedesededesovesedesedevedededeveds veveseveve 
t 

CALL EXCMS ('CLRSCRN' ) 
WRITE(6, 1008) 

1 
WRITE(6,*) ' READING LAMA AND UGVEX MATRICIES 
RIIEBCE) ey) 


THIS SECTION READS THE LAMA VECTOR AND THE UGVEX 
MATRIX AND STORES THEM IN MEMORY FOR FURTHER RECALL OF 
DESIRED LOCATION DATA. 


READ(4,1001) NAM 
READ(4, 1002)(LAMA(1I) ,1=1,100) 
READ(4,1001) NAM 
DO 5 J = 1,100 

READ(4, 1002) (UGVEX(I,J) , I=1, 684) 
CONTINUE 


FORMAT( 1X, A6) 
FORMAT( 1X,8E15. 8) 
FORMAT( 1X, ////) 


CALL EXCMS ('CLRSCRN' ) 

Seveskeseslcvlestesicvevesviese STARTING MODE NUMBER PY rise ivelori wr iseyorkseisrlarisc arise) 
*% SMODE 7 TO 100 (INTEGER) “**% 

SMODE=10 


WRITE (16,700) SMODE 


FORMAT (' ','STARTING MODE NUMBER: ',1I2) 


Vesesovevesesleedevevesevevese Sevevesesevevestevedevevede 


NUMBER OF MODES TO SCAN 
*#* MODE 1 TO 93 (INTEGER) 


MODE= 3 
EMODE = SMODE + MODE - 1 


WRITE (16,701) MODE 


FORMAT (' ','NUMBER OF MODES SCANNED: ',12) 


Sevesevevesesevesededeve NOISE INPUT POSITION devodetesedtedevedevededtedede sci se 
*v NODE 1 TO 114 (INTEGER) (IF 0 THEN NO NOISE INPUT) 


NODE= 8 
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GMA01080 
GMA0 1090 
GMA01100 
GMA01110 
GMA01120 
GMA01130 
GMA01140 
GMA01150 
GMA01160 
GMA01170 
GMA01180 
GMA01190 
GMA01200 
GMA01210 
GMA01220 
GMA01230 
GMA01240 
GMA01250 
GMA0 1260 
GMA01270 
GMA01280 
GMA01290 
GMA01300 
GMA01310 
GMA01320 
GMA01330 
GMA01340 
GMA01350 
GMA01360 
GMA01370 
GMA01380 
GMA01390 
GMA01400 
GMA01410 
GMA01420 
GMA01430 
GMA01440 
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GMA01570 
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mOz 


C2 0) Ca 


© 


op J app sl > Dal ip 


GaGa Co 


WRITE (16,702) NODE 
FORMAT (' ',' NOISE NODE LOCATION: el) 


MWaesedevetedsedededevede SAMPLING TIME KItetedevetevessteveds sXe 
“we SAMPT MUST BE LESS THAN OR EQUAL TO SAMPTM ** 
SAMPT = .05 
SAMPTM = ((2. ODO*PI)/SQRT(LAMA(EMODE) ))/2. ODO 
IF (SAMPT. GE. SAMPTM) THEN 
SAMPT=SAMPTM 
ENDIF 


WRITE (16,900) MIN 
FORMAT (' ‘,2X,'MIN: ',F8. 3) 
WRITE (16,703) SAMPT 


FORMAT (' ','SAMPLING TIME: ',D12.4) 


eri Roses eel alse) ae aig) y DAMPING FACTOR rior aeise ik aclar ok ieee ae 
vv DAMP 0.0 TO 1.0 CREAL*8) 
DAMP=. 01 


WRITE (16,704) DAMP 


FORMAT (' ','DAMPING FACTOR: ',D12.4) 


weer PLANT NOISE VARIANCE **% 
we PNVARX, PNVARY, PNVARZ GT 0.0 


SF1=2. 5D06 


PNVARX=1. ODOO*SF1 
PNVARY=1. ODOO*SF1 
PNVARZ=1. ODOO*SF1 


veveve MEASUREMENT NOISE VARIANCE *** 
v%* MNVARX, MNVARY, MNVARZ GT 0.0 
SF=1. 0 

MNVX1=5. 5D-03*SF 

MNVY1=5. 5D-03*SF 

MNVZ1=5. 5D-03*SF 


CALL EXCMS ('CLRSCRN' ) 
WRITE (6,1008) 


WRITE (6,*) ' PROGRAM RUNNING' 
Diss ernr COO Aki aes NOISE INPUT LOCATION Sevetedevevotetedededde 
ROWN3 = NODE*6 

ROWN2 = (NODE*6) - 1 

ROWN1 = (NODE*6) - 2 

COUNT = 0 
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45 
40 


58 
60 


MIaanwqn 


aQan 


90 


2 
94 


rites geile ayeiselae ses 


INITIALIZE MATRICIES 
DO 40 I =1,3 
DO 45 J = 1,3 
RMN(I,J)=0. 0 
CONTINUE 
CONTINUE 


DO 47 I=1,50 
DO 46 J=1,50 
IDENT(1,J)=0. 0 
PR( 1 (J )=0. 0 
CONTINUE 
CONTINUE 


DO 48 K=1,50 
IDENT(K,K)=1. 0 
CONTINUE 


week TNITIALIZE RMN AND QPN MATRICES *%*%* 


Deseo I=1,3 
DO 58 J=1,3 
QPN(I,J)=0. 0 
CONTINUE 
CONTINUE 


RMN( 1, 1)=MNVX1*"**2 
RMN( 2, 2)=MNVY 17"*2 
RMN(3,3)=MNVZ1#*2 
QPN( 1, 1)=PNVARX***2. 0 
QPN( 2, 2)=PNVARY*"*"2. 0 
QPN( 3, 3)=PNVARZ**2. 0 


FORMAT (' ',’ ") 


aloalen!s slevledlesledleslovlese vleslevevedletliedevededevevevedededededesevetevedevey Sfodstentecianiaew oul walantentsulou! Severe secies! alsulwulon’ 
ves vier is FV EV FLED EL ER GR ER EL EDEL ER EL EL. ET Go AT FR GR ELEN EV Ed.EV_Ersd. Ed EBs OX river fs ever és €% 7d eaves @#ye @ e @e even ex ey de ad rivirty eve. ri 


Saas aa 


Soule ute ale a!on' sale n'snlselsulen'sulauteule seule a's nls uleulswaeulenlsalenleules's uses nloulsnloulesnlontsulsu'sutsntealon! Se verleclosle=) we wloule uta ule uteutontoalenlou'e 
CV EVEN EVEN EN ELEVEN GL ED EVEN EN ENED OVEN ELEVEN EC ERED EN EVEL EVEL ENED EN FLED ER EVEN EVEN EVER PRES CD EV ERED EDEN EV CL EVEN EV UNEVEN ED 


dere resed BEGIN MAIN PROGRAM 


qe 


CALL STMTRX(EMODE , SMODE , SAMPT, DAMP, PHI,GAMMA,A,B,LAMA,UGVEX,C, 
+ ROWN1,ROWN2,ROWN3, BN) 


were PRE-LOOP PORTION OF KALMAN FILTER 
JK=SMODE*2-2 


M2=2**MODE 
DO 94 I=1,3 
DO 92 J=1,M2 
JL=IK+J 
SUM=0. 0 
DO 90 K=1,3 
SUM=SUM+QPN(1,K)*BN( JL, K) 
CONTINUE 
TMP1(J,1)=SUM 
CONTINUE 
CONTINUE 
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GMA02330 
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a7 
98 


100 


C) C) 2 


Joa 


DO 98 I=1,M2 
IL=JIK+I 
DO 97 J=1,M2 
SUM=0. 0 
DO 96 K=1,3 
SUM=SUM+BN( JL,K)*TMP1(J,K) 
CONTINUE 
BQBT(I,J)=SUM 
CONTINUE 
CONTINUE 


M2=2*MODE 

DO 100 I=1,M2 
DO 99 J=1,M2 
TMP3(1I,J)=0. 0 
CONTINUE 

CONTINUE 

JL=JIK+M2 

DO 9375 I=1,3 
DO 9374 J=1,JL 
C(1,J)=C(1,J)*SF 
CONTINUE 

CONTINUE 


Sedesesevesdedecdedesetdestceseskevesededeseseslcdeveskedesdcdesdcakvdvdkvoededkdcskvakcakcclesdke sede desevede sede veskeseskeseskeoesesk 
wa nlealeniea's 

UIE THIS SECTION COMPUTES THE STATE UPDATE 
sevevesetesovevesesescvevevovevedcvestesededededvevedededededevevededcvdevesvesedvedcvededesevedetevevededededetedcdevede 


ESUM=0. 0 


s‘e et alsatsa'e 
SeVev at a8 


SETS LOOP FOR THE ITERATIONS NECESSARY TO OBSERVE = *%ve%% 
THE SYSTEM ePORK THE-NUNBER OF MINUIES SEEGIE TED 


ri eri wi eri erty 


rier iwi eri wrry 


LOOP = INT(C(MIN*60. 0)/SAMPT) 
PRT=(DBLE( LOOP) )/1200. 0 
PRTA=(DBLE( LOOP) )/2400. 0 
CNTA=0. 0 

ACHK=1. 0D-10 

H1=0. 
H2=0. 
H3=0. 
H4=0. 
H5=0. 
H6=0. 
TCHK=MIN*60. 0 
CONTINUE 


Oo @O.O'oQ 0 © 


TIME = TIME+ SAMPT 


CGN=CGN+1. 0 
CNTA=CNTA+1. 0 
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170 


180 


LSS 
190 


C2) CIC) 


wiek START OF KALMAN FILTER *%% 
M2=2*MODE 
weve COMPUTATION OF PK*AT %*%% 


JK=2*SMODE -2 
DO 175 I=1,M2 
DO 170 J=1,M2 
JL=IK+J 
SUM=0. 0 
DO 165 K=1,M2 
JM=JK+K 
SUM =SUM+PK(I,K)*A(JL, JM) 
CONTINUE 
TMP3(1I,J)=SUM 
CONTINUE 
CONTINUE 


wie’ COMPUTATION OF A*(PK*AT)+ BQBT = PK1 *%%* 


DO 190 I=1,M2 
JL=JK+I 
DO 185 J=1,M2 
SUM=0. 0 
DO 180 K=1,M2 
JM=JK+K 
SUM=SUM+A( JL, JM)*TMP3(K, J) 
CONTINUE 
PK1(1,J)=SUM+BQBT(I,J) 
CONTINUE 
CONTINUE 


Sev aeververereve 


“wee COMPUTE PK1*CT wvsrr% 


DO 205 I=1,M2 
DO 200 J=1,3 
SUM=0. 0 
DO 195 K=1,M2 
JM=JK+K 
SUM=SUM+PK1(1,K)*C(J,JM) 
CONTINUE 
TMP1(1I,J)=SUM 
CONTINUE 
CONTINUE 


Tevedevevedesedeslovevede se 


veteve COMPUTE C**( PK1“CT )+RMN Teds ve 
DO 220 I=1,3 
DO 215 J=1,3 
SUM=0. 0 
DO 210 K=1,M2 
JM=JK+K 
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210 


zal 
220 


202 €- crc2 C2 ©) 


Zoo 


240 
245 


el Il aaa 


SUM=SUM+C( I ,JM)*TMP1(K,J) 


CONTINUE 
TMP2(1,J)=SUM+RMN( I,J) 
CONTINUE 

CONTINUE 


wx COMPUTATION OF THE INVERSE OF C*PK1*CT + R 


CALL DLINRG ( 3,7MP2,3,TMP2,3) 


wee COMPUTE CT*INV(C*PK1*CT+R) 


DO 245 I=1,M2 
SL=IK+I 
DO 240 J=1,3 
SUM=0. 0 
DO 235 K=1,3 


SUM=SUM+C(K, JL)*TMP2(K,J) 


CONTINUE 
TMP1(1I,J)=SUM 
CONTINUE 

CONTINUE 


CRP Ce eT ee Se 
ee ie tro bite hic ok bk torres 


vie COMPUTE PK1*C*¥ INV(C*PK1*CT+R ) = G rr%etedek 


DO 260 I=1,M2 
DO 255 J=1,3 
SUM=0. 0 
DO 250 K=1,M2 


SUM=SUM+PK1(1,K)*TMP1(K,J) 


CONTINUE 
G(1,J)=SUM 
CONTINUE 

CONTINUE 


N9=DABS((G(1,1)-H1)/G(1,1)) 
IF (N9. GT. ACHK) THEN 
Gore 7377 

END IF 
N9=DABS((G(1,3)-H2)/G(1,3)) 
IF (N9. GT. ACHK) THEN 

Goro 7377 

END IF 
N9=DABS((G(2,1)-H3)/G(2,1)) 
IF (N9. GT. ACHK) THEN 
GO TO 7377 

END IF 
N9=DABS((G(2,3)-H4)/G(2,3)) 
IF (N9. GT. ACHK) N 
Goro 7377 

END IF 
N9=DABS((G(3,3)-H5)/G(3,3)) 
IF (N9. GT. ACHK) THEN 
GOTO 7377 
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Be Ta 


265 


270 
275 
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END IF 

N9=DABS((G(M2 ,3)-H6)/G(M2,3)) 
IF (N9. GT. ACHK) THEN 

GO TO 277 

ENDTIE 

GO TO 400 


CONTINUE 

Hi=COied) 
H2=G(1,3) 
H3=G(2,1) 
H4=G(2,3) 
H5=G( 3,3) 
H6=G(M2,3) 


PROC ick’ LE PINE) Sane 
GO TO 400 

END IF 

Ti SCCGN GEESE) THEN 


WRITE (6,*) 'TIME= ', TIME, ' SEC. ' 


WRITE (6,*) 'N9= ', NO 
CGN=0. 0 
END 


vette COMPUTE IDENT - G*C 


DO 275 I=1,M2 
DO 270 J=1,M2 
JL=IK+J 
SUM=0. 0 
DO 265 K=1,3 
SUM=SUM+G( I ,K)*C(K, JL) 
CONTINUE 
TMP3(1,J)= IDENT(1,J)-SUM 
CONTINUE 
CONTINUE 


wf. = ='s a eit ee - a 
Seseseveveveveslevedesevesevedevetevedede 


weve COMPUTE PK= (INDENT - G*C)*PK1 


DO 290 I=1,M2 

DO 285 J=1,M2 

SUM=0. 0 
DO 280 K=1,M2 
SUM=SUM+TMP3( I ,K)*PK1(K,J) 
CONTINUE 

PKC 1I,J)=SUM 

CONTINUE 


CONTINUE 
Tedededsvedetedededededededededesese dese 
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C 
400 
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7154 
e155 


CJ-C2 C2 C26 


C2 G2" O) 


END OF KALMAN FILTER PART 1 - START OF PART 2 veer 


GO TO 9991 
CONTINUE 


WRITE (20,1008) 

WRITE (20,*) 'TIME= ',TIME 

DO 384 I=1,M2 

Wee 2055350) GC1,1),G(1,2),G(I,3) 
CONTINUE 

HORHAiG | 5X,D15.8 ,5X,D15.8 ,5X,D15.8 ) 
WRITE (20,*) 'N9= ',NQ 


werk COMPUTE AGC = A - G*C 


M2=2*MODE 
JK=2*SMODE -2 


DO 7155 I=1,M2 
JL=JK+I 
DO 7154 J=1,M2 
JIM=JIK+J 
SUM=0. 0 
DO 7153 K=1,3 
SUM=SUM+G( I ,K)*C(K,JM) 
CONTINUE 
AGC(I,J)=A( JL, JM) -SUM 
CONTINUE 
CONTINUE 


we COMPUTE THE EIGENVALUES OF AGC 

CALL DEVCRG (M2, AGC, 100, EVAL, EVEC, 100) 
woe PRINT EVAL (EIGENVALUE) MATRIX 

DO 7157 I=1,M2 


WRITE (20,*) ‘I= ', I, 'EIG= ', EVAL(I) 
CONTINUE 


STOP 
END 


Fesccedededevesedevesesvescucsesesesesevede skeak seve destcdedevededevedeve sede se ves sede seve seve tose seve se seslese seve ak ese ve ve ok ae 
THIS SUBROUTINE COMPUTES THE STATE TRANSITION MATRIX FOR EACH 
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Carga Ct) C) Co 


ee 1 o> 1 pm J a | 


He 


OF THE 100 MODES 


SevedevesedctetetetetedoleseleleidteiviciinchkniiccccclvhRehekicckhkeTchiddedcdcKeikded 


SUBROUTINE STMTRX(EMODE , SMODE ,T,D, PHI ,GAMMA,A,B,LAMA,UGVEX,C, 
ROWN1 , ROWN2 , ROWN3, BN) 


REAL*8 WN,GMA,PHI(2,2,100),GAMMA( 2,100) ,EGT,T,COSW1T,SINW1T 
REAL#87 W1,D,AC 200, 200) | B@2Z0072 ete. 20007 BN Coe! 

REAL LAMA( 100), UGVEX( 684, 100) 

INTEGER SMODE,R,EMODE,JJ,KK,ROWN1,ROWN2, ROWN3 


1) COO 6S a , 100 
WN = DBLE(SQRT(LAMA(I))) 
GMA = D*WN/2. 0 
EGT = DEXP( -GMA*T) 
W1 = DSQRT( (WN*2) -(GMA%*2) ) 
COSW1T = DCOS(W1*T) 
SINW1T = DSIN(W1*T) 


IF(WN. EQ. 0) THEN 


PHI(1,1,1) = EGT*COSW1T 

PHI(1,2,1) =T 

PHI(2,1,1) = 0 

PHI(2,2,1) = EGT*COSW1T 

GAMMA(1,1I) = 0 

GAMMA(2,1I) = 0 
ELSE 
PHI(1,1,1) = EGI*(COSW1T + (GMA*(W1**(-1)))*SINW1T) 
PHI(1,2,1) = (W1*(-1))*EGT*SINW1T 
PHI(2,1,1) = -(WN**2)*(W1**(-1))*EGI*SINW1T 
PHI(2,2,1) = EGT*(COSW1T - (GMA*(W1**(-1)))*SINWIT) 


GAMMA( 1, I)=CWN****( -2) )*( 1. DO-EGT*COSW1T-EGT*(GMA/W1)*SINW1T) 
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gqdilaqni © Cleo ClO) Gc) qaan 


qaqagIaQqanrIiangdany 


C20) ©) C9 02 0) CG) Gs ore 


GAMMA(2,1) = (W1**(-1))*EGT*SINW1T 


ENDIF 


CONTINUE 


R= 1 


DO 610 K=1 , 100 


A(R,R) = PHI(1,1,K) 
A(R,R+1) = PHI(1,2,K) 
AC(R+1,R) = PHI(2,1,K) 
A(R+1,R+1) = PHI(2,2,K) 


“=< B MATRIX FOR MULTIPLYING CONTROL TORQUES 


B(R,1) = GAMMA(1,K)*DBLE(UGVEX(412,K)) 
B(R,2) = GAMMA(1,K)*DBLE( UGVEX(413,K)) 
B(R,3) = GAMMA(1,K)*DBLE(UGVEX(414 ,K)) 
B(R+1,1) = GAMMA(2,K)*DBLE(UGVEX(412,K)) 
B(R+1,2) = GAMMA(2,K)*DBLE(UGVEX(413,K)) 
B(R+1,3) = GAMMA(2,K)*DBLE(UGVEX(414,K)) 


wee BN MATRIX FOR MULTIPLYING THE NOISE DISTURBANCES 


BN(R,1)=GAMMA( 1,K)*DBLE( UGVEX(ROWN1 ,K)) 
BN(R,2)=GAMMA( 1 ,K)**DBLE( UGVEX(ROWN2,K)) 
BN(R,3)=GAMMA( 1,K)**DBLE( UGVEX( ROWN3 ,K) ) 
BN(R+1, 1)=GAMMA( 2, K)*DBLE(UGVEX(ROWN1 ,K) ) 
BN(R+1,2)=GAMMA( 2 ,K)**DBLE( UGVEX( ROWN2 ,K) ) 


47 


GMA06070 
GMA06080 
GMA06090 
GMA06100 
GMA06110 
GMA06120 
GMA06130 
GMA06140 
GMA06150 
GMA06160 
GMA06170 
GMA06180 
GMA06190 
GMA06200 
GMA06210 
GMA06220 
GMA06230 
GMA06240 
GMA06250 
GMA06260 
GMA06270 
GMA06280 
GMA06290 
GMA06300 
GMA06310 
GMA06320 
GMA06330 
GMA06340 
GMA06350 
GMA06360 
GMA06370 
GMNA06380 
GMA06390 
GMA06400 
GMA06410 
GMA06420 
GMA06430 
GMA06440 
GMA06450 
GMA06460 
GMA06470 
GMA06480 
GMA06490 
GMA06500 
GMA06510 
GMA06520 
GMA06530 
GMA06540 
GMA06550 
GMA06560 
GMA06570 
GMA06580 
GMA06590 
GMA06600 
GMA06610 
GMA06620 


C2 C5 CGC Cl fa C2 Ga Gao Co) 2 ae a 


DE Se 


PApca Cl Ox) 


BN(R+1,3)=GAMMA( 2, K)*DBLE( UGVEX( ROWN3, K) ) 


R = R+2 
CONTINUE 


wee C MATRIX PRODUCTION *%**%* 


jJJ=-1 
DO 640 I=1,100 
JJ=II+1 
KK=1+JJ 


C(1,KK) 
CCZeRK) 
C(3,KK) 


DBLE( UGVEX(418 ,1)) 
DBLE( UGVEX( 419 ,1)) 
DBLE( UGVEX(420,I)) 


KK=KK+1 

C(1,KK)=0. 0 
C(2,KK)=0. 0 
C(3,KK)=0. 0 


CONTINUE 


RETURN 
END 
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GMA06630 
GMA06640 
GMA06650 
GMA06660 
GMA06670 
GMA06680 
GMA06690 
GMA06700 
GMA06710 
GMA06720 
GMA06730 
GMA06740 
GMA06750 
GMA06760 
GMA06770 
GMA06780 
GMA06790 
GMA06800 
GMA06810 
GMA06820 
GMA06830 
GMA06840 
GMA06850 
GMA06860 
GMA06870 
GMA06880 
GMA06890 
GMA06900 
GMA06910 
GMA06920 
GMA06930 
GMA06940 
GMA06950 
GMA06960 
GMA06970 
GMA06980 
GMA06990 
GMA07000 
GMA07010 
GMA07020 
GMA07030 
GMA07040 
GMA07050 
GMA07060 


Co) Caer? ClO) Co Cr GC? O20) 1) 0 


oo «- 


Careared CC) C2 CIC Care 


APPENDIX  B. KALMAN OBSERVER AND PLANT SIMULATION 


Sentetewavevewavtetevewanions ra onto = 
BETES TE TE TE ETE TE PETE TE BCT UCR T ETS ToT E POPE ETE TELE I VESTED TEER IT ELECT ESN V CT SEPSIS 


Pevevedeve S IMRUN Wevewese re 
‘eeeevert =~ ADAPTED TO READ KALMAN FILETER G MATRICE Ca a 
wavest THEN RUN ALL N- MODES OF THE PLANT WHILE soi ae 
‘ieee ~USING A KALMAN FILTER TO OBSERVE M Edd 
*kkAX NUMBER OF STATES een 


WIRE Ree RaoakekkKecsdlecceecectsetest 


UTR TRIER REE RRR RIC RRR Raa RRR om RRE REE, RR; 
esede rete VARIABLE DECLARATIONS Week 


RRR RRR RRR Ron ok eelvclcelelelclclelelelelelelese 


EXTERNAL STMTRX,EXCMS 
CHARACTER*6 NAM 
CHARACTER*1 AGAIN,CORECT,RAGAIN 


INTEGER ROWN1,ROWN2 ,ROWN3, COUNT, NODE ,MODE ,KQ,EMODE, SMODE ,R2M,C2M 
INTEGER CT,CF,KADJ,CFADJ,LOOP, PRNT,JJ,JK,N1,JR,KR,MR, ISEED ,M2 


INTEGER ITYPE(200), IPVT(100), NS, NF, SN, FN 
INTEGER JL,J1,JM , JP, JQ, KA, KB, KC, KD, KE, KF, KG 


REAL LAMA(100), UGVEX(684, 100) ,RNODE,RMODE ,MIN 

REAL*8 PHI(2,2,100),GAMMA( 2,100) ,EGT,GMA,WN,W1,X1T,X2T,TIME 
REAL*8 A(200,200),B(200,3),F(3, 50), IMPLSE,ENERGY 

REAL*8 COSW1T,SINW1T,X(200) 

REAL*8 TCX,TCY,TCZ,DAMP,SAMPT, PI ,SAMPTM,SUM1,SUM2,SUM3,SUMC 
REAL*8 C(3,200), IDENT( 50, 50), RMN(3,3), QPN(3,3) 

REAL*8 Y(3) , BN(200,3) 

REAL*8 PNVARX, PNVARY, PNVARZ 

REAL*8 MNVARX, MNVARY, MNVARZ 

REAL*8 SUM, RNDM(6), RND1, RND2 

REAL*8 XH( 50) ,BOBT( 50, 50) 

REAL*8 SF1 

REAL**8 TMP1( 50,3), TMP2(3,3), TMP3( 50, 50) 

REAL*8 G( 50,3) 

REAL*8 XH1( 50) ,DY(3) , ES,ED,ESUM,CGN, PRT 

REAL*8 WT , WD(3), BNWD(200) 

REAMeeex( 200) , V03), SF ueTO, CIT, ESS 


Rete ose to, SDE eh2(100),.4DELL, ERS, PREI, E3SC100), xXS( 100) 


Seveveslevedevevoteseve tes sede se dete seve dete te teve Fesece FRR RIT RRR RT RRR CIC 


auscacscs VARIABLE DEFINITIONS Weew es 
Pere dede He Hee We Pew Weed He HHH AN KIRIN NIA IN INIA INLINE 


STMTRX = SUBROUTINE EXTABLISHES STATE TRANSITION MATRICIES 
LAMA = VECTOR OF THE SQUARE OF THE NATURAL FREQUENCIES 
UGVEX = MODE POSITONS AND SLOPES OF THE NODAL POINTS 

PHI = STATE TRANSITION MATRICIES FOR EACH MODE 
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SIM00010 
SIMO0020 
SIM00030 
SIM00040 
SIMO00050 
SIM00060 
5IM00070 
5 IMO00080 
SIMO0090 
5IM00100 
5SIM00110 
SIMO00120 
5IM00130 
5IM00140 
5SIM00150 
SIM00160 
»IM00170 
SIM00180 
SIM00190 
SIMO0200 
SIM00210 
SIM00220 
5 IM00230 
SIM00240 
SIM00250 
SIM00260 
SIM00270 
SIM00280 
SIMO00290 
SIM00300 
SIM00310 
SIM00320 
SIM00330 
SIM00340 
SIM00350 
SIM00360 
SIM00370 
SIM00380 
SIM00390 
SIM00400 
SIM00410 
SIM00420 
SIM00430 
SIM00440 
SIM00450 
SIM00460 
SIM00470 
SIM00480 
SIM00490 
SIM00500 
SIMO00510 


CG GC ea Cape GG) Co Ca Cl) Ga CD C2 Cy C2 CaeG2 C20 a .Co GG Gs Ca > CCG) Ga aS) CG Ga) C2). Cor) Co) CGC) Ce Cantera) C) CCD Ca C2 


GAMMA = INPUT TRANSITION MATRIX 

A = DIAGONAL MATRIX CONSISTING OF PHI 

B = INPUT MATRIX OF GAMMA AND CONTROL SLOPES 
DAMP = DAMPING FACTOR 

SAMPT = SAMPLING TIME 

TCX, TCY, TCZ = CONTROL TORQUE VALUES 

ENERGY = TOTAL SYSTEM ENERGY 

IMPLSE = IMPULSE INPUT FUNCTION 

MIN = NUMBER OF MINUTES SYSTEM WILL BE OBSERVED 
SMODE = NUMBER OF STARTING MODE CINT) 

MODE = NUMBER OF MODES (INT) 

EMODE = NUMBER OF THE LAST MODE (INT) 

NODE = NUMBER OF THE NOISE INPUT MODE (INT) 
weve NOISE SLOPE LOCATIONS IN DATA MATRIX *** 
ROWN1 = X-SLOPE LOCATION 

ROWN2 = Y-SLOPE LOCATION 

ROWN3 = Z-SLOPE LOCATION 

C = OUTPUT MATRIX FOR Y 

IDENT = IDENTITY MATRIX 

RMN = MEASUREMENT NOISE COVARIANCE MATRIX 
QPN = PLANT NOISE COVARIANCE MATRIX 


PNVARX = PLANT NOISE X-SLOPE VARIANCE 
PNVARY = PLANT NOISE Y-SLOPE VARIANCE 
PNVARZ = PLANT NOISE Z-SLOPE VARIANCE 
MNVARX = MEASUREMENT NOISE X-SLOPE VARIANCE 
MNVARY = MEASUREMENT NOISE Y-SLOPE VARIANCE 
MNVARZ = MEASUREMENT NOISE Z-SLOPE VARIANCE 


ISEED = INITIALIZATION FOR RANDOM NUMBER GENERATOR 

XKAL = X MATRIX 

Y = OUTPUT MATRIX 

RNDM = RANDOM NUMBERS USED FOR WHITE NOISE IN MEASUREMENTS AND 
IN PLANT FORCES 

BN = B MATRIX TO MULTIPLY NOISE DISTURBANCES 

TNX,TNY,TNZ= NOISE TORQUES X,Y,Z SLOPES 

M2=2**MODE 


sealssesestcskesteslesteske sevesicscslokslksicuicsicctcsleat slvatkeslevic ale 


SAMPLE OF SPACE EXEC FILE 


THIS FILE MUST BEGIN IN COLUMN 1 AND RUN WITH THE FOLLOWING 
SEQUENCE FOR THE INITIAL RUN OF THE PROGRAM: 


FORTVS SPACE 
SPACE 
LOAD SPACE (START 


(COMPILES PROGRAM) 
(EXECUTES EXEC FILE) 
(LOADS AND EXECUTES PROGRAM) 


SUBSEQUENT PROGRAM RUNS CAN ELIMINATE "FORTVS SPACE" IF NO 
CHANGES HAVE BEEN MADE TO THE PROGRAM, AND CAN ELIMINATE 
RUNNING THE EXEC FILE. 


FI 4 DISK THESIS INPUT B (PERM 

FI 8 DISK UTILITY DATA (RECEM Ys BLOCK 1 33erukn 

FI 11 DISK CNTRL OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 
FI 13 DISK GAMMA OUTPUT (RECFM VS BLOCK 133 PERM 

FI 14 DISK MODE OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 


SO 


SIM00520 
SIMO00530 
SIMO00540 
SIM00550 
SIMO00560 
SIM00570 
SIM00580 
SIMO00590 


) 


; 


SIM00600 


SIM00610 
SIM00620 
SIM00630 
SIMO00640 
SIM00650 
SIM00660 
SIM00670 
SIM00680 
SIM00690 
SIM00700 
SIM00710 
SIM00720 
SIM00730 
SIM00740 
SIM00750 
SIM00760 
SIM00770 
SIM00780 
SIM00790 
SIMO0800 
SIMO00810 
SIM00820 
SIM00830 
SIMO00840 
SIM00850 
SIM00860 
SIM00870 
SIM00880 
SIM00890 
SIM00900 
SIM00910 
SIM00920 
SIM00930 
SIM00940 
SIM00950 
SIM00960 
SIM00970 
SIM00980 
SIM00990 
SIM01000 
SIM01010 
SIM01020 
SIMO1030 
SIMO01040 
SIM01050 
SIM01060 
SIM01070 


Cry G2 C2 CP CIID 


LD i Se Jl Sip I app Jl | C70) 


CB 1a Se lal > He Be 


FI 16 DISK COST OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 
FI 17 DISK PRT OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 
FI 18 DISK ERROR DATA (RECFM F BLOCK 80 LRECL 80 PERM 
RislgeOrokeeND FILE (RECFM F BLOCK 80 LRECL 80 PERM 
FI 20 DISK GMAT FILE (RECFM F BLOCK 80 LRECL 80 PERM 


devededeveuvuicdvas deve dssvsedesededededededesesdeseskdesedevedededsevesedesededecldesesevevesesevevetededededvese seve vevencreveve 


PARAMETER (JR=5243, KR=5397, MR=262139) 


MIN =1.00 


WT=1. ODOO 
PI = 4.0D0 * ATAN(1. ODO) 


Fee TER Tee HH HK TATA TE I I IIH IAAT EH ICT NTR TET TIENEN VIET TE TEETER TE IE TE TEN TE Te TEE IE 
cae READ LAMA AND UGVEX MATRICIES 


STEM Te Tee Tec Te Te eee Teas Te Fs FoF eve eae Te Te TE TS ETE Te Te TSE TER TEE TE TENE TEE TE Te TEE TEE TE TE TE TERETE TE TEE HE GS 


Sevevevese 


CALL EXCMS ('CLRSCRN' ) 
WRITE(6,1008) 
WRITE(6,*) ' 
WRITE(6,*) ° ' 

THIS SECTION READS THE LAMA VECTOR AND THE UGVEX 
MATRIX AND STORES THEM IN MEMORY FOR FURTHER RECALL OF 
DESIRED LOCATION DATA. 


READ(4,1001) NAM 

READ(4, 1002)(LAMA(I), I=1, 100) 

READ(4,1001) NAM 

Doms! = 1100 
READ(4,1002)(UGVEX(1,J),1=1,684) 

CONTINUE 


FORMAT( 1X, A6) 
FORMAT(1X,8E15. 8) 
Hema (1xX5////) 


CALL EXCMS ('CLRSCRN' ) 


Citrisrisrir rit oiets STARTING MODE NUMBER Jedsdesedsdetetcdedededededs 
%% SMODE 7 TO 100 (INTEGER) **7%%* 


SMODE= 7 


WRITE (16,700) SMODE 


FORMAT (' ','STARTING MODE NUMBER: ‘' ,12) 


SesevetevodedesedsseteKisdede NUMBER OF MODES TO SCAN sevcesesesotslodcvevededve 


** MODE 1 TO 93 (INTEGER) 
MODE=20 


EMODE = SMODE + MODE - 1 


| 


READING LAMA AND UGVEX MATRICIES' 


SIM01080 
SIM01090 
51M01100 
SIMO1110 
SIM01120 
SIMO1130 
SIM01140 
SIM01150 
SIMO1160 
SO 170 
SIM01180 
SIM01190 
SIM01200 
SIM01210 
SIM01220 
SIM01230 
SIM01240 
SIM01250 
SIMO1260 
SIMO1270 
SIM01280 
SIMO1290 
SIM01300 
SIM01310 
SiMG1T320 
SIM01330 
SIM01340 
SIM01350 
SIM01360 
SIM01370 
SIM01380 
Sieg 159'0 
SIM01400 
5IM01410 
SIM01420 
SIM01430 
SIM01440 
SIM01450 
SIM01460 
SIM01470 
SIM01480 
SIM01490 
SIM01500 
Simos 10 
el ylehls 210) 
STMO1>30 
SIM01540 
SIM01550 
SIM01560 
SIM01570 
SIMO1580 
Send 20 
SIM01600 
SIM01610 
SIM01620 


Bo9 


[ep] 


C2) C2 C2 C2 C2 


WRITE (16,701) MODE 


FORMAT (' ','NUMBER OF MODES SCANNED: ‘',12) 


r Oiischiche iwi sc forint Gets NOISE INPUT POSITION crac hd isriarierigcis cise aioe oe erisel 
*%* NODE 1 TO 114 CINTEGER) (IF O THEN NO NOISE INPUT) 


NODE= 8 


WRITE (16,702) NODE 


FORMAT (' ',' NOISE NODE LOCATION: ESD) 
wei START AND STOP FOR PLANT 

SN=7 

FN=20 

NS=SN* 2-1 

NF=SN*2+2*FN-2 

WRITE (16,899) SN,FN 

FORMAT (' ','PLANT -- SN= ',I5,' FN= ',I5) 


Sevlevedeslevevevededesevove 


SAMPLING TIME dededededevedededededededeve 
wie SAMPT MUST BE LESS THAN OR EQUAL TO SAMPTM ** 
SAMPT = 0.05 
SAMPTM = ((2. ODO*PI)/SQRT( LAMA(CEMODE)))/1.0D01 
IF (SAMPT. GE. SAMPTM) THEN 
SAMPT=SAMPTM 
ENDIF 


WRITE (16,900) MIN 
FORMAT Ce 2s, MEN@eee roe 
WRITE (16,703) SAMPT, SAMPTM 


FORMAT (' ','SAMPLING TIME: ',D12.4,2X,'SAMPTM= ',D15. 8) 


8 L) toe! 
vevleveslosedeveslevetledtevesedese Sesesieveveseveseseveveseseds 


DAMPING FACTOR 
7 DAMP O70) TOs). OF CREAT 3) 
DAMP=. 01 


WRITE (16,704) DAMP 


FORMAT (' ','DAMPING FACTOR: ',D12.4) 


‘te PLANT NOISE VARIANCE ** 

*% PNVARX, PNVARY, PNVARZ GT 0.0 
SF1=2. 5D06 

SF=1. ODOO 


PNVARX=1. ODOO*SF1 
PNVARY=1. ODOO*SF1 
PNVARZ=1. ODOO*SF1 


weve MEASUREMENT NOISE VARIANCE *** 
*%* MNVARX, MNVARY, MNVARZ GT 0.0 
MNVARX=1. OD-03 *SF 

MNVARY=1. OD-03 *SF 


Sy 


SIM01630 
SIM01640 


SIM01650. 
SIMO01660_ 


SIM01670 
SIM01680 
SIMO1690 
SIMO01700 


SIMO1710' 


SIMO17 2¢ 
S ENO? a 


SIMO1740 


SIM01750 
SIMO1760 
SIMO1770 
SIM01780 
SIMO01790 
SIM01800 
SIM01810 
SIM01820 
SIM01830 
SIM01840 
SIMO01850 
SIM0O1860 
SIMO1870 
SIM01880 
SIM01890 
SIM01900 
SIMO i1a 0 
SIM01920 
SIM01930 
SIMO1940 
SIMO1950 
SIMO01960 
SIMO1970 
SIM01980 
SIM01990 
SIM02000 
SIMO2010 
SIMO02020 
SIM02030 
SIM02040 
SIM02050 
SIM02060 
SIM02070 
SIMO2080 
SIM02090 
SIMO2100 
SIM02110 
SIM02120 
SIMO02130 
SIM02140 
SEMGZISC 
SIM02160 


{ 


C27 © 


Sp ae LP 1a | 


MNVARZ=1. OD-03 *SF 


WRIte (16,711) 


FORMAT(' ','PLANT NOISE VARIANCE: 


Weis (16,712) 


FORMAT(' ' ,6X,'PNVARX' ,13X,'PNVARY' ,13X,'PNVARZ') 
WRITE (16,713) PNVARX, PNVARY, PNVARZ 
BOmAme 92n,b15.8,2X,H15.8,2X,E15.8) 


WRITE( 16,714) 


FORMAT(' ', ‘MEASUREMENT NOISE: ') 


Wit 6, 715)) 


FORMAT(' ',6X,'MNVARX' ,13X,'MNVARY' ,13X,'MNVARZ') 
WRITE( 16,713) MNVARX,MNVARY , MNVARZ 


CALL EXCMS ('CLRSCRN' ) 
WRITE (6,1008) 
WRITE (6,*) ' 


Verssevevevevevevevess 


ROWN3 
ROWN2 
ROWN1 
COUNT 


NODE**6 
(NODE*6) - 1 
(NODE*6) - 2 
0 


Sevevevevovovevevesieve 
DO 48 K=1,50 


IDENT(K,K)=1. 0 
CONTINUE 


2 
nee. = 0.0 
CONTINUE 


WRITE(C6, 1008) 


WRITE (6,%*) ' INITIALIZE RMN AND QPN MATRICES 
we INITIALIZE RMN AND QPN MATRICES *** 


CONTINUE 


RMN( 1, 1 )=MNVARX*2"2 
RMN( 2, 2)=MNVARY**2 
RMN( 3, 3)=MNVARZ**2 
QPN(1,1)=PNVARX***2 
QPN(2,2)=PNVARY**2 


NOISE INPUT LOCATION 


INITIALIZE MATRICIES 


oe 


PROGRAM RUNNING' 


Pededetevedededede cede 


tesevededesedevededesevedcvededs 


SIMO02170 
SIMO2180 
Sime2T 0 
SIM02200 
SIM02210 
SIM02220 
SIM02230 
SIM02240 
SIM02250 
SIM02260 
SIM02270 
SIM02280 
S1H022310 
SIM02300 
SIM02310 
SIM02320 
SIM02330 
SIM02340 
SIMO02350 
SIM02360 
SIM02370 
SIM02380 
SIM02390 
SIM02400 
SIM02410 
SIM02420 
SIM02430 
SIM02440 
SIM02450 
SIM02460 
SIM02470 
SIM02480 
SIMO02490 
SIM02500 
SIM02510 
SIMO2520 
SIM02530 
SIM02540 
SIM02550 
SIM02560 
SIM02570 
SIM02580 
SIM02590 
SIM02600 
SIM02610 
SIM02620 
SIM02630 
SIM02640 
SIM02650 
SIMO02660 
SIM02670 
SIM02680 
SIM02690 
PEN 0Z7 00 


C2) C2 0) ©) ©) C26) 


ao 


384 


So? [is > Je Ge a ep a 


9777 1 


QPN(3,3)=PNVARZ** 2. 0 


WRITE(6, 1008) 
WRITE(6,*) ‘ ENTER STMTRX 


Sevevetevededesescvevescdevesesedevesovesevevededevevetevevevesevesedete vevevevesvedessse seve seve seve se veveseve seve vere 


vevedededs BEGIN MAIN PROGRAM 


Yeveveveve 


aeveveveverevese ted dedededededededstevedeve seve seus ved tek te texs te te tevevese tested sede Tete Tete ve TF Fe Ke HK 


CALL STMTRX(EMODE , SMODE ,SAMPT , DAMP, PHI ,GAMMA,A,B, LAMA, UGVEX,C, 
+  ROWN1,ROWN2 ,ROWN3, BN) 


WRITE (16,1008) 

DO 11000 I=1,200 
DO 10900 J=1,3 
C(J,I)= C( J,1)*SF 
CONTINUE 

CONTINUE 


vevevs PRE-LOOP PORTION OF KALMAN FILTER 


M2=2*MODE 
JP=2*SMODE-1 
JQ=2**EMODE 
DO 90 I=1,50 
XH(1)=0. 0 
CONTINUE 


DOES 71 1=1. NZ 
READ ((20,**)°GG1, in eGGl2) Cel.) 
CONTINUE 


WRITE (14,1008) 
DO 384 I=1,M2 


WRITE (14,5350)  G(I,1),G(1I,2),G(I,3) 


CONTINUE 


FORMAT (' ',2X,D15.8,2X,D15. 8,2X,D15. 8) 


Sesetevetssededevesekdvedevededededovesededevedeseveveseseslededededesedevevedevevesesevedeseseetvacdedededeseve 


Feveweev THIS SECTION COMPUTES THE STATE UPDATE 


Sevevevese 


vite ibibitintrintininininivinininitrinintetetEi teat ktrikteipinaioieihintateixtEininretkreietei EEE ieee ee 


DO 9771 I=1,100 
E2(1)=0.0 
E3(1)=0. 0 
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SIMO2710 
SIM02720 
SIM02730 
SIM02740 
SIMO02750 
SIM02760 
SIM02770 
SIM02780 
SIM02790 
SIMO2800 
SIM02810 
SIM02820 
SIM02830 
SIM02840 
SIM02850 
SIM02860 
5SIM02870 
SIMO2880 
SIM02890 
SIM02900 
SIM02910 
SIM02920 
5IM02930 
SIM02940 
SIMOZ95¢ 
SIM02960 
SIM02970 
SIMO02980 
SIMO02990 
SIMO03000 
SIM03010 
SIM03020 
SIM03030 
SIMO03040 
SIMO03050 
SIM03060 
SIM03070 
5SIM03080 
SIM03090 
SIM03100 
> IM@aiaG 
SIMOs2o 
SIM03130 
SIM03140 
> INOS i150 
SIM03160 
SIM03170 
5IM03180 
SIM03190 
SIM03200 
SiMO2 2G 
SIM03220 
SIM03230 
SIM03240 
SIM03250 


OQ kr 
S 
= 


Ga Cy Ca tae Ca ca G2 


oo 


CTG=0. 0 


weer =~ SETS LOOP FOR THE ITERATIONS NECESSARY TO OBSERVE 
weewerere =6TTHE SYSTEM FOR THE NUMBER OF MINUTES SPECIFIED 
WRITE (6,1008) 

WRITE (6,*) ' START STATE UPDATE 

LOOP = INT((MIN*60. 0)/SAMPT) 

PRT= (DBLE( LOOP) )/30. 0 

Boo 


sevesedese 
Sertedeveve 


DO 400 L = 0, LOOP 
TIME = DBLE(L)*SAMPT 


IF(L. EQ. 0) THEN 

IMPLSE =(1. OD06*SF1)/(DSQRT( SAMPT) ) 
ELSE 

IMPLSE = 0. 0DO 
ENDIF 


TO=0. 0 
‘erkvevesk, RANDOM NUMBER GENERATOR deteecest 


DO 101 I=1,6 

ISEED=MOD( ISEED*\JR+KR , MR) 
RND1=(DBLE( ISEED)+0. 5D00) /DBLE( MR) 

ISEED=MOD( ISEED**JR+KR , MR) 
RND2=(DBLE( ISEED)+0. 5D00) /DBLE(MR) 

RNDM( 1 )=DSQRT( -2. 0*DLOG(RND1) )*DCOS(6. 2831853D00*RND2) 
CONTINUE ; 


sevevevevevesesevedesedesevedeseseaicss 


Grr=CTi+1, 0 
seveckevs START OF STATE UPDATE *** 
* X°200 


Pac OMPUEn AX°200 = A°200 X 200 


“tet COMPUTE AX = AX 


JK=SMODE*2-2 
JP=JK+1 
JQ=2**EMODE 


DO 5015 I=NS,NF 
SUM=0. 0 
DO 5010 K=NS,NF 
SUM=SUM+A(1,K)*X(K) 
CONTINUE 
AX(1)=SUM 
CONTINUE 


we COMPUTE WD°3 
WD( 1)=PNVARX*RNDM( 1)*TO+IMPLSE 


WD( 2)=PNVARY*RNDM( 2)*TO 
WD(3)=PNVARZ*RNDM( 3) *TO 


5D 


SIM03260 
SIM03270 
Siis2260 
SIM03290 
SIM03300 
SIM03310 
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G2) C1 


C2 G2 C1" Coe) 
) 
{ 
=) 


mca 


5045 


050 


wee COMPUTE BNWD°200 =BN°200 X 3. * WD°3 


DO 5035 I=NS,NF 
SUM=0. 0 
DO 5030 K=1,3 
SUM=SUM+BN(I,K)*WDCK) 
CONTINUE 
BNWD( 1)=SUM 
CONTINUE 


wee" COMPUTE X°200 =AX°200 + BNWD°200 


DO 5040 I=NS,NF 
X(I)= AX(I) + BNWD(T) 
IF (DABS(X(I)).LT. 1.0D-60) THEN 
X(I)=1. OD-60 
END IF 


CONTINUE 
9% COMPUTE V°3 
V(1)=MNVARX**RNDM(4) 


VC 2)=MNVARY*RNDM(5 ) 
V(3)=MNVARZ**RNDM(6 ) 


vee COMPUTE Y°3 = C°3 X 200 * xX°200 + V°3 


DO 5050 I=1,3 
SUM=0. 0 
DO 5045 K=NS,NF 
SUM=SUM+C( I ,K)*X(K) 
CONTINUE 
Y(I)=SUM+V(I) 
CONTINUE 


eorlorleovlevlesloclovoveclonte Jedeuloadavoue 
TATA rd var ty es #8 #8 Oe ee rd vir or # i #¥ es FN eR ed ed €% 


weve START OF KALMAN FILTER **** 


M2=2*MODE 


week COMPUTE XH1 = A*XH 


DO 300 I=JP,J0 

SUM=0. 0 
DO 295 K=JP,JQ 
SUM=SUM+A(I,K) * XH(K) 
CONTINUE 

XH1(1)=SUM 

CONTINUE 


Serle eels ale Seclodledlevioudewocioua 
Ae av ae ev a8 ee ON aN aS at Al 68 ae a8 
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CaG779 


310 


1s 


IOIAAAW 


320 


ho 
U1 


CC? 2 C7) C907 Cd WD OD 


340 


2100 


*#% COMPUTE DY = Y - C*¥XH1 


DO 315 I=1,3 

SUM=0. 0 
DO 310 K=JP,JQ 
SUM=SUM+C(I ,K)*XH1(K) 
CONTINUE 

DY(I)=Y(1)-SUM 

CONTINUE 


Vededededodedetedevededtetodedetere 


were COMPUTE XH = XH1 + G*DY 


DO 325 I=1,M2 

J1=JP-1+I 

SUM=0. 0 
DO 320 K=1,3 
SUM=SUM+G( I ,K)*DY(K) 
CONTINUE 

XH( J1)=XH1(J1)+SUM 

IF (DABS(XH(J1)). LT. 1. 0D-60) THEN 
XH(J1)=1. 0*D-60 

END IF 


CONTINUE 


dessdeders END OF KALMAN ROUTINES  wceecedreess 


wee COMPUTATION OF ESUM **%* 
DO 340 I=JP,JQ 

XDEL= X(1I)-XH(1) 
XDEL1=XDEL*XDEL*SAMPT 
E2(1)=E2(1)+XDEL1 
XSC1I)=XSCI)+XC1)*XC1I)*SAMPT 
E3( 1I)=E2(1)/XS(1) 

CONTINUE 


CGN=CGN+1. 0 
DE Gewine® 1. 0. OR. CGN. GIT. PRT) THEN 


WRITE (6,*) ‘TIME= ', TIME, ' SEC. ' 


WRITE (17,1008) 
WRITE (16,1008) 
WRITE (16,2100) TIME 


WRITE (17,2100) TIME 

FORMAT(' ','TIME= ',F9. 3) 

DO 380 I=JP , JQ 

WRITE (16,4500) I,X(I) ,I ,XH(I) 
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SIM04380 
SIM04390 
SIM04400 
SIM04410 
SIM04420 
SIM04430 
SIM04440 
SIM04450 
SIM04460 
SIM04470 
SIM04480 
SIM04490 
SIM04500 
SIM04510 
SIM04520 
SIM04530 
SIM04540 
SIM04550 
SIM04560 
SIM04570 
SIM04580 
SIM04590 
SIM04600 
SIM04610 
SIM04620 
SIM04630 
SIM04640 
SIM04650 
SIM04660 
SIM04670 
SIM04680 
SIM04690 
SIM04700 
SIM04710 
SIM04720 
SIM04730 
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SIM04750 
SIM04760 
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SIM04780 
SIM04790 
SIM04800 
SIMC4810 
SIM04820 
SIM04830 
SIM04840 
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SIM04860 
SIM04870 
SIM04880 
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SIM04920 


WN 
\O 
\O 


Gi-C). £20) CC) OC) 


Gy 2G) C) 


WE 10h ayo 3 Obl 5 2 ui) 
CONTINUE 


s9E3C(I) =, « XS(T) 


CGN=0. 0 
END IF 
FORMAT (' 
FORMAT (' 


y, XC',13;, J=  yDI5.8,2%, SHO (135 peo 
5%, 150%5 ce loeen 


CONTINUE 


DO 401 I=JP,JQ 
WRITE (19756530) =i ee2( 1 je eott) xsi) 
CONTINUE 


STOP 
END 


SKiveketevevekdevevesdekeseteveteteveveveteveteteteteveteseKye te vesedesetetedetetevedevete seve dete Teseve Fe Tee eT eT N FETE TET 
THIS SUBROUTINE COMPUTES THE STATE TRANSITION MATRIX FOR EACH 
OF THE 100 MODES 


skicveseveucvestvevekseseseoescstevevesovedsestvedesesesonkokvesevosevesteseskveske seven seskscvesevevesesesesk ves kvevevesede 


SUBROUTINE STMTRX(EMODE, SMODE,T,D,PHI,GAMMA,A,B,LAMA,UGVEX,C, 
ROWN1 , ROWN2 , ROWN3 , BN) 


REAL*8 WN,GMA,PHI(2,2,100) ,GAMMA(2,100),EGT,T,COSW1T,SINW1T 
REAL*8 W1,D,A( 200,200) ,B(200,3),C(3, 200) ,BN( 200, 3) 

REAL LAMA( 100) ,UGVEX( 684, 100) 

INTEGER SMODE,R,EMODE,JJ,KK,ROWN1,ROWN2,ROWN3 


DO 600s es =—1 , 100 
WN = DBLE(SQRT( LAMA(I))) 
GMA = D*WN/2.0 
EGT = DEXP( -GMA*T) 


W1 = DSQRT((WN***2) -( GMA**2) ) 
COSW1T = DCOS(W1*T) 
SINW1T = DSIN(W1*T) 


IF(WN. EQ. 0) THEN 


PHI(1,1,1) = EGT*COSW1T 
PHI(1,2,1) =T 
PHI(2,1,1) = 0 
PHI(2,2,1) = EGT*COSW1T 
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S1IM05310 
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SIM05470 
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gqoaaand 


qolaalaqand >) IaAANRIMdHd aqanan mqaqan gqaan 


Iaqaaiadknna 


GAMMA(1,I) = 0 

GAMMA(2,1) = 0 
ELSE 
PHI(1,1,1) = EGT*(COSW1T + (GMA%*(W1*( -1)))*SINW1T) 
PHI(1,2,I1) = (W1%*(-1))*EGT*SINWIT 
PHI(2,1,1) = -(WN%*2)%*(W1%*(-1))*EGT*SINWIT 
PHI(2,2,1) = EGT*(COSW1T - (GMA*(W1*(-1)))*SINWIT) 


GAMMA( 1, I)=CWN***( -2))**(1. ODO-EGT*COSW1T -EGT*(GMA/W1)*SINW1T) 
GAMMA(2,I) = (W1**(-1))*EGT*SINW1T 


ENDIF 


CONTINUE 


= il 


DO 610 K=1 , 100 


ACR,R) = PHI(1,1,K) 
A(R, R+1) = PHI(1,2,K) 
A(R+1,R) = PHI(2,1,K) 
A(R+1,R+1) = PHI(2,2,K) 


wee B MATRIX FOR MULTIPLYING CONTROL TORQUES 


B(R,1) = GAMMA(1,K)*DBLE(UGVEX( 412 ,K)) 
B(R,2) = GAMMA(1,K)*DBLE(UGVEX( 413 ,K)) 
B(R,3) = GAMMA(1,K)*DBLE(UGVEX(414,K)) 


B(R+1,1) = GAMMA(2,K)*DBLE(UGVEX(412 ,K)) 


a2 
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C2 C2 CG GaGa) CC Ca 


maaanina 


a" 
© 


Le I > 1 I Ie a 1 9 


De Rl > 1 > 


BCR+1,2) 


GAMMA( 2 ,K)*DBLE( UGVEX( 413 ,K) 
BiG lero) K) 


GAMMA( 2 ,K)*DBLE(C UGVEX( 414, 


a ed 


wie BN MATRIX FOR MULTIPLYING THE NOISE DISTURBANCES 


BN(CR, 1)=GAMMA(1,K)**DBLE( UGVEX(CROWN1,K) ) 
BN(R,2)=GAMMA(1,K)*DBLE( UGVEX( ROWN2,K) ) 
BN(CR,3)=GAMMA(1,K)*DBLE( UGVEX( ROWN3 , K) ) 
BN( R+1,1)=GAMMA( 2,K)*DBLE( UGVEX( ROWN1,K) ) 
BN(R+1, 2)=GAMMA( 2 ,K)*DBLE( UGVEX( ROWN2, K) ) 
BNC R+1,3)=GAMMA( 2 ,K)*DBLE( UGVEX( ROWN3 , K) ) 


R = R+2 
CONTINUE 


vee C MATRIX PRODUCTION *** 


JJ=-1 

DO 640 I=1,100 

JJ=II+1 

KK=I el 

C(1,KK) = DBLE(UGVEX(418,1I)) 

C(2,KK) = DBLE(UGVEX(419,I)) 

C(3,KK) = DBLE(UGVEX(420,1)) 
=KK+1 

C(1,KK)=0. 0 

C(2,KK)=0. 0 

C(3,KK)=0. 0 

CONTINUE 
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APPENDIX C. PROGRAM TO ESTIMATE NOISE IN KALMAN FILTER 


C2 c2 C2 C2 Co C2 Ca Ca CD Cd Cc 


ap Il Gy Fs ap 


Ci€2 Ca Ca GC) Cl C.-C) 2 OC) 


FROM UNOBSERVED MODES 


dededededededcdedededededededesededs sede dededk ede de vede Fete aR etete Re RRERREERERER RII 
Seaesesteve SPAC 24 tevedaedk 


‘eeeee ~~ ADAPTED TO RUN N MODES OF THE PLANT AND aedecedesk 
‘eevee = COMPUTE THE NOISE IN THE KALMAN FILTER serknledese 


‘eevee = FROM THE UNOBSERVED MODES serkcedere 
Sesesedtkescdevedsacdededeslesesesevesesesesescaededesesdededesetedesedededesededk Teves se deste reve Feet RK 


Seakesesesesdestcsededeseucdededesesesesesesesdedecdcsedesesedevdedesede descdeveddescaedesesetededkeacsedededesedededesesedk 


HRA VARIABLE DECLARATIONS HaeH ev 
sestcstesestcvleslescdeacsesesestcscslevcveslcsicsevesedesededcsededescdedesesedededededesesesesededcwdesdcsdedesesevesdeseseskesede 


EXTERNAL STMTRX ,EXCMS 

CHARACTER*6 NAM 

CHARACTER*1 AGAIN, CORECT,RAGAIN 

INTEGER ROWN1,ROWN2, ROWN3,COUNT, NODE ,MODE ,KQ,EMODE , SMODE ,R2M,C2M 
INTEGER CT,CF,KADJ,CFADJ,LOOP,PRNT,JJ,JK,N1,JR,KR,MR, ISEED,M2 
INTEGER NO, NS, NF, SN, FN 

INTEGER JL,J31,5M , SP, JG. KA, KB KG. KO KE NE KG 


REAL LAMA( 100), UGVEX(684,100),RNODE,RMODE ,MIN 
REAL*8 PHI(2,2,100),GAMMA(2,100),EGT,GMA,WN,W1,X1T,X2T, TIME 
REAL*8 A(200,200),B(200,3),F(3, 50), IMPLSE,ENERGY 
REAL*8 COSW1T,SINW1T,X(200) 

REAL*8 DAMP, SAMPT, PI ,SAMPTM,SUM1,SUM2,SUM3 ,SUMC 
REAL*8 C(9,200), RMN(3,3), QPN(3,3) 

REAL*8 BN(200,3) 

REAL*8 PNVARX, PNVARY, PNVARZ 

REAL*8 MNVARX, MNVARY, MNVARZ 

REAL*8 SUM, RNDM(6), RND1, RND2 

REAL*8 ES,ED,ESUM,CGN, PRT 

REAL*8 WT , WD(3), BNWD(200), EX1(9) 

REAL*8 EX(9),AX(200) , SF , TO, CTT, ESS 

REAL*8 CTG, XDEL, XDEL1, ERS, PRT1 

REAL*8 SF1 


KivdekkkkdckkeaekkkdedkdksededededksdescdedcsededesckdedcslekkkkKcKRKR RRR, 


ena VARIABLE DEFINITIONS SEE IOSs 
HIRAM AAA HAMAR IANA KIA IKI INN NIT ANNA IARI ITN IRIAN I IRI IT 


STMTRX = SUBROUTINE EXTABLISHES STATE TRANSITION MATRICIES 
LAMA = VECTOR OF THE SQUARE OF THE NATURAL FREQUENCIES 
UGVEX = MODE POSITONS AND SLOPES OF THE NODAL POINTS 

PHI = STATE TRANSITION MATRICIES FOR EACH MODE 

GAMMA = INPUT TRANSITION MATRIX 

A = DIAGONAL MATRIX CONSISTING OF PHI 
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SPA00010 
SPA00020 
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SPA00080 
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Cieea G2 G2 C2 C2 C2 C2 CCI CD Crerrerera aa? 0) I) 0) CF Or ara 0 a rar I ar a Os Oe a Ss Oe Ga Ss a Gs Os a a) 


B = INPUT MATRIX OF GAMMA AND CONTROL SLOPES 
DAMP = DAMPING FACTOR 

SAMPT = SAMPLING TIME 

IMPLSE = IMPULSE INPUT FUNCTION 

MIN = NUMBER OF MINUTES SYSTEM WILL BE OBSERVED 
SMODE = NUMBER OF STARTING MODE (INT) 

MODE = NUMBER OF MODES (INT) 

EMODE = NUMBER OF THE LAST MODE (INT) 

NODE = NUMBER OF THE NOISE INPUT MODE (INT) 
weve NOISE SLOPE LOCATIONS IN DATA MATRIX *** 
ROWN1 = X-SLOPE LOCATION 

ROWN2 = Y-SLOPE LOCATION 

ROWN3 = Z-SLOPE LOCATION 

GC = OUTPUT MATRIX FOR Y 

IDENT = IDENTITY MATRIX 

RMN = MEASUREMENT NOISE COVARIANCE MATRIX 
QPN = PLANT NOISE COVARIANCE MATRIX 


PNVARX = PLANT NOISE X-SLOPE VARIANCE 

PNVARY = PLANT NOISE Y-SLOPE VARIANCE 

PNVARZ = PLANT NOISE Z-SLOPE VARIANCE 

MNVARX = MEASUREMENT NOISE X-=SLOPE VARIANCE 
MNVARY = MEASUREMENT NOISE Y-SLOPE VARIANCE 
MNVARZ = MEASUREMENT NOISE Z-SLOPE VARIANCE - 


ISEED = INITIALIZATION FOR RANDOM NUMBER GENERATOR 

RNDM = RANDOM NUMBERS USED FOR WHITE NOISE IN MEASUREMENTS AND 
IN PLANT FORCES 

BN = B MATRIX TO MULTIPLY NOISE DISTURBANCES 


FeterereWr IH ere I SAMPLE OF SPACE EXEC FILE 
THIS FILE MUST BEGIN IN COLUMN 1 AND RUN WITH THE FOLLOWING 
SEQUENCE FOR THE INITIAL RUN OF THE PROGRAM: 


FORTVS SPACE 
SPACE 
LOAD SPACE (START 


(COMPILES PROGRAM) 
(EXECUTES EXEC FILE) 
(LOADS AND EXECUTES PROGRAM) 


SUBSEQUENT PROGRAM RUNS CAN ELIMINATE "“FORTVS SPACE" IF NO 
CHANGES HAVE BEEN MADE TO THE PROGRAM, AND CAN ELIMINATE 
RUNNING THE EXEC FILE. 


FI 4 DISK THESIS INPUT (PERM 

FI 8 DISK UTILITY DATA (RECFM VS BLOCK 133 PERM 

FI 11 DISK CNTRL OUTPUT (CRECFM F BLOCK 80 LRECL 80 PERM 
FI 13 DISK GAMMA OUTPUT (RECFM VS BLOCK 133 PERM 

FI 14 DISK MODE OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 

FI 16 DISK COST OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 

FI 17 DISK PRT OUTPUT (RECFM F BLOCK 80 LRECL 80 PERM 

FI 18 DISK ERROR DATA (RECFM F BLOCK 80 LRECL 80 PERM 

FI 19 DISK END FILE (RECFM F BLOCK 80 LRECL 80 PERM 

FI ZO DISK GMAT FILE CRECFM F BLOCK 80 LRECL 80 PERM 


seicdevescsdesestedeseteseakvcdestesededeskeskacvdesetesesc ted toskevlede teseskieseskeseskeakskskscskeskseskskakscskesksesbotskscakstsestets 
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Seseseskoesetesedesesedevesesesesesec 


SPA00500 
SPA00510 
SPA00520 
SPA00530 
SPA00540 
SPA00550 
SPA00560 
SPA00570 
SPA00580 
SPA00590 
SPA00600 
SPA00610 
SPA00620 
SPA00630 
SPA00640 
SPA00650 
SPA00660 
SPA00670 
SPA00680 
SPA00690 
SPA00700 
SPA00710 
SPA00720 
SPA00730 
SPA00740 
SPA00750 
SPA00760 
SPA00770 
SPA00780 
SPA00790 
SPA00800 
SPA00810 
SPA00820 
SPA00830 
SPA00840 
SPA00850 
SPA00860 
SPA00870 
SPA00880 
SPA00890 
SPA00900 
SPA00910 
SPA00920 
SPA00930 
SPA00940 
SPA00950 
SPA00960 
SPA00970 
SPA00980 
SPA00990 
SPA01000 
SPA01010 
SPA01020 
SPA01030 
SPA01040 
SPA01050 


C2 C2: C2 Cl 


202) @ 


702 


PARAMETER (JR=5243, KR=5397, MR=262139) 


MIN 15 30 


WT=1. ODOO 
PI = 4.0D0 * ATAN(1. 0D0) 


Vere Fore He Fe vee re Vee ve Fe Foye Fe Fe Fe He TENE He NE He He He FEN eH TTC TI TIN HIT IK IEEE RENE RE 1 


AIH READ LAMA AND UGVEX MATRICIES HARI 
Seseceveseslesecesteueskeslese dese deseskeslesesestcleste ves alent clkesesteskvokslentatcleslosescseslestcc dese seats ake de vesentedescsesesk 


CALL EXCMS ('CLRSCRN' ) 
WRITE(6, 1008) 
WRITE(6,*) ' 
WRITE(6,*) ' ' 

THIS SECTION READS THE LAMA VECTOR AND THE UGVEX 
MATRIX AND STORES THEM IN MEMORY FOR FURTHER RECALL OF 
DESIRED LOCATION DATA. 


READ(4,1001) NAM 
READ(4,1002)( LAMA(I) , I=1, 100) 
READ(4,1001) NAM 
DO 5 J = 1,100 

READ(4, 1002) (UGVEX(1I,J),I=1,684) 
CONTINUE 


FORMAT( 1X, A6) 
FORMAT( 1X, 8E15. 8) 
FORMAT(1X,////) 


CALL EXCMS ('CLRSCRN' ) 


oriaciscise nitrite STARTING MODE NUMBER sevevesevevetetesedesevesere 
%% SMODE 7 TO 100 CINTEGER) **** 


SMODE= 17 


WRITE (16,700) SMODE 
FORMAT (' ','STARTING MODE NUMBER: ‘',I2) 
Jevededededsdededededsdese 


Heseskccavceokskesbaevdesededese NUMBER OF MODES TO SCAN 
** MODE 1 TO 93 (INTEGER) 


MODE=3 
EMODE = SMODE + MODE - 1 


WRITE (16,701) MODE 
FORMAT (' ','NUMBER OF MODES SCANNED: ‘',12) 


ederedere deeded ye NOISE INPUT POSITION HAA A RAKAKAKII 


*% NODE 1 TO 114 (INTEGER) (IF 0 THEN NO NOISE INPUT) 
NODE= 8 


WRITE (167702 @ NODE 


FORMAT (' ',' NOISE NODE LOCATION: ' ,I5) 
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READING LAMA AND UGVEX MATRICIES' 


SPA01060 

SPA01070 
SPA01080 
SPA01090 
SPA01100 © 
SPA01110 | 
SPA01120 — 
SPA0113G 
SPA01140 
SPA01150 
SPA01160 
SPA01170 
SPA01180 
SPA01190 
SPA01200 
SPA01210 
SPA01220 
SPA01230 
SPA01240 
SPA01250 
SPA01260 
SPA01270 
SPA01280 
SPA01290 
SPA01300 
SPA01310 
SPA01320 
SPA01330 
SPA01340 
SPA01350 
SPA01360 
SPA01370 
SPA01380 
SPA01390 
SPA01400 
SPA01410 
SPA01420 
SPA01430 
SPA01440 
SPA01450 
SPA01460 
SPA01470 
SPA01480 
SPA01490 
SPA01500 
SPA01510 
SPA01520 
SPA01530 
SPA01540 
SPA01550 
SPA01560 
SPA01570 
SPA01580 
SPA01590 
SPA01600 
SPA01610 
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714 


weivdedkewoeseacdedeke START AND STOP FOR PLANT 
SN=17 
FN=4 
NS=SN*2-1 
NF=SN*2+2*FN-2 
WRITE (16,899) SN,FN 
FORMAT (' ','PLANT -- SN= ',15,' FN= ',I5) 
Flvkedevedevede sede ded SAMPLING TIME Hasek Hed Kx 
w#%* SAMPT MUST BE LESS THAN OR EQUAL TO SAMPTM ** 
SAMPT = 0.05 
SAMPTM = ((2. ODO*PI)/SQRT( LAMA(EMODE)))/1.0D01 
IF (SAMPT. GE. SAMPTM) THEN 
SAMPT=SAMPTM 
ENDIF 


WRITE (16,900) MIN 
FORMAT (' ',2X,'MIN: ',F8.3) 
WRITE (16,703) SAMPT, SAMPTM 


FORMAT (' ','SAMPLING TIME: ',D12.4,2X,'SAMPTM= ',D15. 8) 


Sedesvedkdedededkdeseucdesedte DAMPING FACTOR Kistedetedededederksedeverk 
*#* DAMP 0.0 TO 1.0 CREAL*8) 
DAMP=. 01 


WRITE (16,704) DAMP 

FORMAT (' ','DAMPING FACTOR: ',D12.4) 
NO=3 

veveck PLANT NOISE VARIANCE **%* 

** PNVARX, PNVARY, PNVARZ GT 0.0 


SF=1. ODO 

SF1=2. 5D06 
PNVARX=1. ODOO*SF1 
PNVARY=1. ODOO*SF1 
PNVARZ=1. ODOO*SF1 


wave MEASUREMENT NOISE VARIANCE **% 
weve MNVARX, MNVARY, MNVARZ GT 0.0 
MNVARX=1. 0D-03 *SF 

MNVARY=1. OD-03 *SF 

MNVARZ=1. OD-03 *SF 


WRITE (16,711) 

FORMAT(' ','PLANT NOISE VARIANCE: ') 

WRITE (16,712) 

FORMAT(' ',6X,'PNVARX' ,13X,'PNVARY' ,13X, 'PNVARZ') 
WRITE (16,713) PNVARX, PNVARY, PNVARZ 

BOmiame " .2X%,E15.8,2X,E15. 8y2X3E15. 8) 

WRITE( 16,714) 

FORMAT(' ','MEASUREMENT NOISE: ° ) 
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SPA01620 
SPA01630 
SPA01640 
SPA01650 
SPA01660 
SPA01670 
SPA01680 
SPA01690 
SPA01700 
SPA01710 
SPA01720 
SPA01730 
SPA01740 
SPA01750 
SPA01760 
SPA01770 
SPA01780 
SPA01790 
SPA01800 
SPA01810 
SPA01820 
SPA01830 
SPA01840 
SPA01850 
SPA01860 
SPA01870 
SPA01880 
SPA01890 
SPA01900 
SPA01910 
SPA01920 
SEAL ZO 
SPA01940 
SPA01950 
SPA01960 
SPA01970 
SPA01980 
SEBO 1990 
SPA02000 
SPA02010 
SPA02020 
SPA02030 
SPA02040 
SPA02050 
SPA02060 
SPA02070 
SPA02080 
SPA02090 
SPA02100 
SPA02110 
SPA02120 
SPA02130 
SPA02140 
SPA02150 
SPA02160 
SPA02170 


Co C2 


Sp I ap Be a | 


58 
60 


CC) G)..C) <a 


D> a 


10900 
11000 
C 


+ 


WRITE( 16,715) 


FORMAT(' ',6X,'MNVARX' ,13X, 'MNVARY' ,13X, 'MNVARZ' ) 
WRITE( 16,713) MNVARX,MNVARY , MNVARZ 


CALL EXGMS ( CLRSCRN » 


WRITE (6,1008) 
WRITE (6,*) ' 


Sesedesevesedesc sexes 


NOISE INPUT LOCATION 


ROWN3 = NODE*6 

ROWN2 = (NODE*6) - 1 
ROWN1 = (NODE*6) - 2 
COUNT = 0 


Taek Hevea sex 


DO 54 K = 1, 200 
X(K) = 0.0 
CONTINUE 


DO 60 I=1,3 
DO 58 J=1,3 
RMN(I,J)=0. 0 
QPN(I,J)=0. 0 
CONTINUE 
CONTINUE 


RMN( 1, 1)=MNVARX***"2. 
RMN( 2, 2)=MNVARY***"2. 
RMN(3, 3)=MNVARZ""2. 
QPN( 1, 1)=PNVARK***2. 
QPN(2,2)=PNVARY"*"2. 
QPN(3,3)=PNVARZ*""2. 


Sevedededededededesesesedesedesesedededededsdedesesesededededededesesesiedesdetedededeenddededetedtondsedededetke ted 


BEGIN MAIN PROGRAM 


Sewsatau'suleuteatsaloutentoatoutou'snts else's alentseloetse'sntoulsntselouisutacionteelsulsulselon'sulantoaleutaaiautseleienisalsetsnion's alaals uals alse), sevesleslesleslsste 
RQ EV FL TL FL FEDERER TEDL FH GR HS ER EL EL FH F% FRE CD COED EX HED riy #u eh 40 ry #2 EO 4% ed FB 4B Fd FD seas #v ee aces #u 48 ees av 4s en ¢ OO 4B Od 4d 4d 


Sedesescse 


CALL STMTRX(CEMODE , SMODE ,SAMPT, DAMP, PHI ,GAMMA,A,B,LAMA,UGVEX,C, 
ROWN1 , ROWN2 , ROWN3, BN) 


WRITE (6,1008) 


WRITE(6,*) ' EXIT STMTRX - - - PRE-LOOP KALMAN’ 


WRITE (6,*) 'COMPUTING C TIMES SF FOR NEW C' 


WRITE (16,1008) 
DO 11000 I=1,200 


DO 10900 J=1,NO 
C(I, 1)= CGD) sr 


CONTINUE 
CONTINUE 


INITIALIZE MATRICIES 


0 
0 
0 
0 
0 
0 


PROGRAM RUNNING' 


Tededededededesesewieg 


TivteaaKeK RK KR 


SPA02180 
SPA02190 
SPA02200 
SPA02210 
SPA02220! 
SPA02230 
SPA02240 | 
SPA02250 
SPA02260 
SPA02270 | 
SPA02280 
SPA02290 
SPA02300 
SPA02310 
SPA02320 
SPA02330 
SPA02340 
SPA02350 

SPA02360 
SPA02370 
SPA02380 
SPA02390 
SPA02400 
SPA02410 
SPA02420 
SPA02430 
SPA02440 
SPA02450 
SPA02460 
SPA02470 
SPA02480 
SPA02490 
SPA02500 
SPA02510 
SPA02520 
SPA02530 
SPA02540 
SPA02550 
SPA02560 
SPA02570 
SPA02580 
SPA02590 
SPA02600 
SPA02610 
SPA02620 
SPA02630 
SPA02640 
SPA02650 
SPA02660 
SPA02670 
SPA02680 
SPA02690 
SPA02700 
SPA02710 
SPA02720 
SPA02730 


Cn Gy-C2 C2 ip al Sp is ap 


pe 1 Ge 
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C2C2 CG) Cy Ce) & 


wre ~PRE-LOOP PORTION OF KALMAN FILTER 


JK=SMODE*2 -2 
M2=2*MODE 


Heaacdtedkeaesdedeveseaeskeskeskakescskvslcstcsedeskeskskcsescslkesesk 


M2=2*MODE 
JP=2*SMODE-1 
JQ=2**EMODE 


DO 8813 I=1,3 
EX(1I)=0. 0 
CONTINUE 


Kaede raver eke seseskseseske desks scsesede ses sesesesedeakc akc akaesvededescskeacscsescvesedke seaweed aevewe weakest 
WIcHIew THIS SECTION COMPUTES THE STATE UPDATE HANI 


acdedevcscdedcveacseseveseucsedevescacdcvesedcdevescdesdvscacakvsdcveveskcvicvevescdedcscslcdeveslesicvcakvcscvcdeslevesesicve 


seveacvdesgs 


ever enene 


SETS LOOP FOR THE ITERATIONS NECESSARY TO OBSERVE 
THE SYSTEM FOR THE NUMBER OF MINUTES SPECIFIED 
(6,1008) 

WRITE (6,%*) ' START STATE UPDATE 

LOOP = INT((MIN*60. 0)/SAMPT) 

PRT= (DBLE(LOOP))/30. 0 

PRT1=(DBLE( LOOP) )/50. 00 

Cik—00 


DO 400 L = 0, LOOP 
TIME = DBLE(L)*SAMPT 


IF(L. EQ. 0) THEN 

IMPLSE =(1. OD06*SF1)/(CDSQRTCSAMPT) ) 
ELSE 

IMPLSE = 0. ODO 
ENE 


TO=0. 0 
wee RANDOM NUMBER GENERATOR ‘ie 


DO 101 I=1,6 
ISEED=MOD( ISEED* JR+KR, MR) 

RND1=( DBLE( ISEED)+0. 5D00) /DBLE( MR) 
ISEED=MOD( ISEED**JR+KR , MR) 
RND2=(DBLE( ISEED)+0. 5D00) /DBLE( MR) 
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SPA02740 
SPA02750 
SPA02760 
SPA02770 
SPA02780 
SPA02790 
SPA02800 
SPA02810 
SPA02820 
SPA028 30 
SPA02840 
SPA02850 
SPA02860 
SPA02870 
SPA02880 
SPA02890 
SPA02900 
SPA02910 
SPA02920 
SPA02930 
SPA02940 
SPA02950 
SPA02960 
SPA02970 
SPA02980 
SPA02990 
SPA03000 
SPA03010 
SPA03020 
SPA03030 
SPA03040 
SPA03050 
SPA03060 
SPA03070 
SPA03080 
SPA03090 
SPA03100 
SPA03110 
SPA03120 
SPA03130 
SPA03140 
SPA03150 
SPA03160 
SPA03170 
SPA03180 
oP A031TI0 
SPA03200 
SPA03210 
SPA03220 
SPA03230 
SPA03240 
SPA03250 
SPA03260 
SPA03270 
SPA03280 
SPA03Z90 


i) 
_ 


ImalIaaagan ‘one 


20 


5010 


5015 


mqoaacq 


qQaa 


© 
{ 
© 


C2 (2 C2 2 Canta cD 


RNDM( I)=DSQRT( -2. 0*DLOG( RND1) )*DCOS( 6. 2831853D00*RND2 ) 


CONTINUE 

WIC ae We ee Tee CII ITE ee 
CTT=CTT+1. 0 

vere START OF STATE UPDATE **%* 


kirk COMPUTE AX°200 = A° 200 X 20Gmecasx 200 


wee COMPUTE AXH = A*XH 


JK=SMODE*2-2 
JP=JK+1 
JQ=2*EMODE 


DO 5015 I=NS,NF 
SUM=0. 0 
DO 5010 K=NS,NF 
SUM=SUM+A( I ,K)**X(K) 
CONTINUE 
AX(1)=SUM 
CONTINUE 


exe COMPUTE WD°3 


WD( 1)=PNVARX**RNDM( 1)**TO+IMPLSE 
WD( 2)=PNVARY*RNDM( 2 )*TO 
WD(3)=PNVARZ*RNDM( 3 )**TO 


*% COMPUTE BNWD°200 =BN°200 X 3. * WD°3 


DO 5035 I=NS,NF 
SUM=0. 0 
DO 5030 K=1,3 
SUM=SUM+BN(I ,K)*WD(K) 
CONTINUE 
BNWD( 1 )=SUM 
CONTINUE 


wes COMPUTE X°200 =AX°200 + BNWD°200 
DO 5040 I=NS,NF 
X(I)= AX(1) + BNWD(I) 
IF (DABS(X(I)). LT. 1.0D-60) THEN 
X(I)=1. OD-60 
END IF 


CONTINUE 


Ler liride iol vestrlaciseisriseise sci oriselseiseisels 


wee START OF KALMAN FILTER *%*% 
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SPA03300 
SPA03310 
SPA03320 
SPA03330 
SPA03340 
SPA03350 
SPA03360 
SPA03370 
SPA03380 
SPA03390 
SPA03400 
SPA03410 
SPA03420 
SPA03430 

SPA03440 ~ 
SPA03450 
SPA03460 
SPA03470 
SPA03480 
SPA03490 
SPA03500 
SPA03510 
SPA03520 
SPA03530 
SPA03540 
SPA03550 
SPA03560 
SPA03570 
SPA03580 
SPA03590 
SPA03600 
SPA03610 
SPA03620 
SPA03630 
SPA03640 
SPA03650 
SPA03660 
SPA03670 
SPA03680 
SPA03690 
SPA03700 
SPA03710 
SPA03720 
SPA03730 
SPA03740 
SPA03750 
SPA03760 
SPA03770 
SPA03780 
SPA03790 
SPA03800 
SPA03810 
SPA03820 
SPA03830 
SPA03840 
SPA03850 


400 


O12) a © 


JK=SMODE* 2-2 
JP=JK+1 
JQ=2*EMODE 
M2=2**MODE 


JL=JQ+1 

DO 8888 I=1,NO 

SUM=0. 0 
DO 8887 K=JL,NF 
SUM=SUM+C(I,K)*X(K) 
CONTINUE 

EX( I )=SUM*SUM*SAMPT+EX( I) 

CONTINUE 


CGN=CGN+1. 0 
eco ba BOs 1.0. ORSCGN. GT. PRT) THEN 


WRITE (16,1008) 
WRITE (16,*) "TIME = ', TIME 


DO 380 I=JP , JQ 

WRITE (16,4500) I,X(I) 

CONTINUE 

FORMAT (' ',2X,'X(',14,')= ',D15.8) 


CGN=0. 0 
END iF 


CONTINUE 

WRITE (11,*) 'SMODE = ', SMODE 
WRITE (11,*) 'EMODE = ', EMODE 
WRITE (11,*) ‘SN = 
WRITE (11,*) 'FN = ',FN 


JL=JQ+1 
DO 9499 I=1,NO 


WRITE -Cll.*) “EX °,Ir ,° xc) 


CONTINUE 


CALL EXCMS ('CLRSCRN’ ) 
WRITE (6,1008) 


STOP 
END 
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SPA03860 
SPA03870 
SPA03880 
SPA03890 
SPA03900 
SPA03910 
SPA03920 
SPA03930 
SPA03940 
SPA03950 
SPA03960 
SPA03970 
SPA03980 
SPA03990 
SPA04000 
SPA04010 
SPA04020 
SPA04030 
SPA04040 
SPA04050 
SPA04060 
SPA04070 
SPA04080 
SPA04090 
SPA04100 
SPA04110 
SPA04120 
SPA04130 
SPA04140 
SPA04150 
SPA04160 
SPA04170 
SPA04180 
SPA04190 
SPA04&200 
SPA04210 
SPA04220 
SPA04230 
SPA04240 
SPA04250 
SPA04260 
SPA04&270 
SPA04280 
SPA04290 
SPA04300 
SPA04310 
SPA04320 
SPA04&330 
SPA04340 
SPA04350 
SPA04&360 
SPA04370 
SPA04380 
SPA04390 
SPA04400 
SPA04410 
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Seatevicstcscvedesededkesesesedesese ake sesestcaeskeskvtaesdeseseseseseskeakeskeskstsesedcdeseseocseslcsesdeakdededesesestesesestsleseakcscsesvese 


THIS SUBROUTINE COMPUTES THE STATE TRANSITION MATRIX FOR EACH 
OF THE 100 MODES 


FevevededevedevedesodedededesetededeseseseketeskeiesetekeeedeiodesedesededeiesededeseseseseseicktelesedesesektedKsesedetove 


SUBROUTINE STMTRX(EMODE,SMODE,T,D,PHI,GAMMA,A,B,LAMA,UGVEX,C, 
i ROWN1, ROWN2 , ROWN3 , BN) 


REAL*8 WN,GMA,PHI(2,2,100) ,GAMMA(2,100) ,EGT,T,COSW1T,SINW1T 
REAL*8 W1,D,A(200,200) ,B(200,3),C(9,200) ,BN(200,3) 

REAL LAMA(100) ,UGVEX(684, 100) 

INTEGER SMODE,R,EMODE,JJ,KK,ROWN1,ROWN2,ROWN3, NN(9), N9, NO 


WRITE (6,*) 'INSIDE STMTRX - - COMPUTE WN, GMA, EFT, W1' 


DO 600 I = 1 ,100 
WN = DBLE(SQRT( LAMA(1))) 
GMA = D*WN/2. 0 
EGT = DEXP( -GMA*T) 
W1 = DSQRT( (WN***2) -(GMA**2)) 
COSW1T = DCOS(W1*T) 
SINW1IT = DSIN(W1*T) 


IFCWN. EQ. 0) THEN 


PHI(1,1,1) = EGT*COSW1T 

PHI(1,2,1) =T 

PHI(2,1,1) = 0 

PHI(2,2,1) = EGT*COSW1T 

GAMMA(1,I) = 0 

GAMMA(2,I) = 0 
ELSE 
PHI(1,1,I1) = EGT*(COSW1T + (GMA*(W1**(-1)))*SINW1T) 
PHI(1,2,1) = (W1**(-1))*EGT*SINW1T 
PHI(2,1,1) = -(WNe#2)*(W1**(-1))*EGT*SINWIT 
PHI(2,2,1) = EGT*(COSW1T - (GMA*(W1**(-1)))*SINWIT) 


GAMMA( 1, 1)=(WN**(-2))*(1. ODO-EGT*COSW1T -EGT*(GMA/W1)*SINW1T) 
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SPA04420 © 


SPA04430 
SPA04440 


SPA04450 | 


SPA04460 
SPA04470 
SPA04480 
SPA04490 
SPA04500 
SPA04510 
SPA04520 
SPA04530 
SPA04540 
SPA04550 
SPA04560 
SPA04570 
SPA04580 
SPA04590 
SPA04600 
SPA04610 
SPA04620 
SPA04630 
SPA04640 
SPA04650 
SPA04660 
SPA04670 
SPA04680 
SPA04690 
SPA04700 
SPA04710 
SPA04720 
SPA04730 
SPA04740 
SPA04750 
SPA04760 
SPA04770 
SPA04780 
SPA04790 
SPA04800 
SPA04810 
SPA04820 
SPA04830 
SPA04840 
SPA04850 
SPA04860 
SPA04870 
SPA04880 
SPA04890 
SPA04900 
SPA04910 
SPA04920 
SPA049 30 
SPA04940 
SPA04950 
SPA04960 





Sp Ws ae som (ol oe ik > ©2 ©2 C27 Ov C2 C2 C2 CacGy Ci-C2 


LP IAT Se Ia 6 aE ep 1 op kale fl gp 


C2 C2) C2 GaGa Ci Ga Co G2) C7 Ci Ca 


GAMMA( 2,1) = (W1**(-1))*EGTI*SINW1T 


ENDIF 


CONTINUE 


WRITE (6,*) 'PHI AND GAMMA COMPUTED' 
WRITE (6,*) ' COMPUTING A, B, BN’ 


R= 1 


DO 610 K = 1 LOC 


ACR,R) = PHI(1,1,K) 
A(R,R+1) = PHI(1,2,K) 
A(R+1,R) = PHI(2,1,K) 
A(R+1,R+1) = PHI(2,2,K) 


“wx B MATRIX FOR MULTIPLYING CONTROL TORQUES 


B(R, 1) 
B(R, 2) 
B(R,3) 
B(R+1,1) 
B(R+1,2) 
B(R+1,3) 


GAMMA( 1,K)*DBLE(UGVEX(412,K)) 
GAMMA( 1,K)*DBLE( UGVEX(413,K)) 
GAMMA(1,K)**DBLE( UGVEX(414,K)) 
GAMMA( 2 ,K)*DBLE(UGVEX(412,K)) 
GAMMA( 2 ,K)*DBLE( UGVEX(413,K)) 
GAMMA( 2, K)*DBLE( UGVEX(414,K)) 


we BN MATRIX FOR MULTIPLYING THE NOISE DISTURBANCES 


BN(R, 1)=GAMMA( 1, K)*DBLE( UGVEX(ROWN1,K)) 
BN(R, 2)=GAMMA( 1,K)**DBLE( UGVEX(ROWN2 ,K)) 
BN(R,3)=GAMMA( 1, K)*DBLE( UGVEX(ROWN3 ,K)) 
BN(R+1, 1)=GAMMA( 2 ,K)**DBLE( UGVEX( ROWN1,K)) 
BN(R+1,2)=GAMMA( 2 ,K)**DBLE( UGVEX(ROWN2 ,K)) 


al 


SPA04970 
SPA04980 
SPA04990 
SPA05000 
SPA05010 
SPA05020 
SPA05030 
SPA05040 
SPA05050 
SPA05060 
SPA05070 
SPA05080 
SPA05090 
SPA05100 
SEe0>110 
SPA05120 
SPA05130 
SPA05140 
SPA05150 
SPA05160 
SPA05170 
SPA05180 
SPA05190 
SPA05200 
SPA05210 
SPA05220 
SPA05230 
SPA05240 
SPA05250 
SPA05260 
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SPA05280 
SPA05290 
SPA05300 
SPA05310 
SPA05320 
SPA05330 
SPA05340 
SPA05350 
SPA05360 
SPA05370 
SPA05380 
SPA05390 
SPA05400 
SPA05410 
SPA05420 
SPA05430 
SPA05440 
SPA05450 
SPA05460 
SPA05470 
SPA05480 
SPA05490 
SPA05500 
SPA05510 
SPA05520 
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BNC R+1,3)=GAMMA( 2 ,K)*DBLE( UGVEX( ROWNS , K) ) 


R = Rt+2 
CONTINUE 


WRITE (6,*) 'A, B, BN COMPUTED' 
WRITE (6,*) 'COMPUTING C' 


week C MATRIX PRODUCTION #¥% 
NO=3 

NN(1)=418 

NN(2)=419 

NN(3)=420 


jJ=-1 

DO 640 I=1,100 

JJ=JJ+1 

DO 9127 K=1,NO 

KK=1+JJ 

NO=NN(K) 

C(K,KK) = DBLE(UGVEX(N9,I)) 
KK=KK+1 

C(K,KK)=0. 0 


CONTINUE 
CONTINUE 


RETURN 
END 
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